{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数据准备"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from matplotlib.pylab import rcParams\n",
    "rcParams['figure.figsize'] = 8, 6\n",
    "pd.options.display.float_format = '{:,.2}'.format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "x = np.array([i * np.pi / 180 for i in range(-180, 60, 5)])\n",
    "y = np.cos(x) + np.random.normal(0, 0.15, len(x))\n",
    "data = pd.DataFrame(np.column_stack([x, y]), columns=['x', 'y'])\n",
    "\n",
    "pow_max = 13\n",
    "\n",
    "# 构造不同幂的x\n",
    "for i in range(2, pow_max):\n",
    "    colname = 'x_%d' % i\n",
    "    data[colname] = data['x']**i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(48, 13)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x119452828>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQaElEQVR4nO3dXYhd13nG8efRyHKhF40qDbGrT4uKENOatDrIgt642G5kUyyS1mDHULupEIGKXuSidRDEJSXgEko/UtFWKCIOiMglpfW0NfgjcfBFO0VnQHH9ETXTgUESbqwoU/cibaXpvL2YrXY8OufMzDnr7K/1/4HQ7LO3z15skmeW3v3utR0RAgC036aqBwAAKAeBDwCZIPABIBMEPgBkgsAHgExsrnoA/Wzfvj327t1b9TAAoFFmZmZ+EBGTvfbVNvD37t2rbrdb9TAAoFFsz/fbR0kHADJB4ANAJgh8AMgEgQ8AmSDwASATBD4AZILABzIzM7+gk6/NamZ+oeqhoGS17cMHkN7M/IKeOD2t64tL2rJ5k84ePaQDe7ZWPSyUhBk+kJHpuWu6vrikpZBuLC5peu5a1UNCiQh8ICOH9m3Tls2bNGHpts2bdGjftqqHhBJR0gEycmDPVp09ekjTc9d0aN82yjmZIfCBzBzYs5WgzxQlHQDIBIEPYCS0eTYHJR0AQ6PNs1mY4QMYGm2ezULgAxgabZ7NQkkHwNBo82wWAh/ASOrc5jkzv8AvoxUIfACtxA3lW1HDB1C6Mlo5uaF8K2b4AEpV1sz75g3lG4tL3FAuJJnh2z5j+z3bb/bZb9t/YnvW9hu2fz7FeQE0T1kz75s3lD/7Sx+hnFNINcP/qqQ/lfS1PvsfkrS/+HOvpD8r/gaQmdQz70E3Zut8Q7kKSQI/Il63vXfAIUckfS0iQtK07Q/ZvjMi3k1xfgDNkbKVkxuzG1NWDX+HpEsrti8Xn30g8G0fk3RMknbv3l3S0ACULdXMu1d5iMDvr1ZdOhFxKiI6EdGZnJysejgAao4nfTemrBn+FUm7VmzvLD4DgKHxpO/GlBX4U5KO2z6n5Zu171O/B+qniU+mcmN2/ZIEvu2vS7pP0nbblyU9I+k2SYqIP5f0oqSHJc1K+pGkX09xXgDp5HQDtIm/2FJI1aXz+Br7Q9JvpjgXgPHI5QZoTr/YVqvVTVsA1cnlBmjOSy6wtAIASfncAM15yQUvV1vqp9PpRLfbrXoYAFqozTV82zMR0em1jxk+gOzk2tlDDR8AMkHgA8AKZazVXxVKOgBQaHvLJjN8ACi0vWWTwAeAQtufRaCkA2BNw7QxNrH1se3PIhD4AAYapq7d5Fp4m1s2KekAGGiYuvY4auFt7p4pCzN8oMEGlU1SlVSGWYpgHO+trfpfDE0sUa1G4AM1spFQGRSCKQNymLp26lp41St51uEXTgoEPlATGw2VQSGYOiCHqWunrIVXveBZ1b9wUiHwgZrYaKgMCsGqAzK1qrtn2nI9WS0TqImbM/ybobLebphx1/CxrKzrOep5Bq2WSeADNUJI5y3FvQKWRwYaos094FjbuO8V0IcPoDZy77Uf99IOzPAB1EJbWh9HMe6b0wQ+0AA51Pbb0vo4qnGW9Qh8oOZymfm2pfWxzgh8oOZymflW3WufAwIfqLmcZr50KY0XgQ/UHDNfpELgAw3AzBcp0IcPAJkg8AFgTOr2IBklHWBMcuidR391bKcl8IExqOP/2VGutdppq5gQEPjAGOTSO4/+BrXTVjUhIPCBMcipdx69DWqnrWpCQOADY0DvfF76lWf6tdNWNSHgBSgAMIJhyzPjquHzAhQAGJNhyzNVPEyXpA/f9mHbF23P2n66x/6nbF+1faH4czTFeQGgauN+aUlKI8/wbU9IOinpQUmXJZ23PRURb6869PmIOD7q+QCgTpp0vyZFSeegpNmImJMk2+ckHZG0OvABoJWastZRipLODkmXVmxfLj5b7Vdsv2H7G7Z39foi28dsd213r169mmBoQP3U7XF75KOstXT+VtLeiLhH0iuSnut1UESciohORHQmJydLGhpQnpsdHX/w8kU9cXqa0EepUgT+FUkrZ+w7i8/+T0Rci4j/LjZPSzqQ4LxA4/Tq6ADKkiLwz0vab/su21skPSZpauUBtu9csfmIpHcSnBdonCZ1dKB9Rr5pGxGLto9LeknShKQzEfGW7S9I6kbElKTfsv2IpEVJP5T01KjnBZqoSR0daB+etAVGwBLIqBuetAXGgCWQ0TS88QqtVEbrIzdg0TTM8NE6Zc28WQIZTUPgo3XKWmucG7BoGgIfrVPmzLspj9QDEoGPFhp25k3HDdqOwEcrbXTmTccNckCXDiA6bpAHAh8QSx4gD5R0ANFxgzwQ+ECBjhu0HSUdAMgEgQ8AmSDwASATBD4AZILAB4BMEPjIThlLJwN1RFsmssISCsgZM3xkZdglFPhXAdqAGT6yMszSyfyrAG1B4CMrwyyhUNYLVYBxI/CRnY0uocCrDNEWBD6wBhZWQ1sQ+MA6sLAa2oAuHQDIBIEPAJkg8AEgEwQ+AGSCwAeATBD4AJAJAh8AMkHgA0AmCHwAyASBDwCZIPABIBNJAt/2YdsXbc/afrrH/tttP1/s/yfbe1OcFwCwfiMHvu0JSSclPSTpbkmP27571WG/IWkhIn5a0h9K+v1RzwsA2JgUM/yDkmYjYi4irks6J+nIqmOOSHqu+Pkbku637QTnBgCsU4rA3yHp0orty8VnPY+JiEVJ70viLRIAUKJa3bS1fcx213b36tWrVQ8HAFolReBfkbRrxfbO4rOex9jeLOknJF1b/UURcSoiOhHRmZycTDA0NMnM/IJOvjarmfmFqocCtFKKN16dl7Tf9l1aDvbHJH1q1TFTkp6U9I+SflXStyIiEpwbLTEzv6AnTk/r+uKStmzepLNHD/GGKSCxkWf4RU3+uKSXJL0j6S8j4i3bX7D9SHHYVyRtsz0r6bOSbmndRN6m567p+uKSlkK6sbik6blb/gEIYERJ3mkbES9KenHVZ59f8fN/SXo0xbnQTof2bdOWzZt0Y3FJt23epEP7uKcPpMZLzFELB/Zs1dmjhzQ9d02H9m2jnAOMAYGP2jiwZytBD4xRrdoyAQDjQ+ADQCYIfDQavfvA+lHDR+3NzC/0vJlL7z6wMQQ+am1QqPfq3Sfwgf4o6aDWBj2QdbN3f8Kidx9YB2b4qLVBD2TRuw9sjOu6pE2n04lut1v1MFAD/Wr4AG5leyYiOr32McNH7fFAFpAGNXwAyASBDwCZIPABIBMEPgBkgsAHgEwQ+ACQCQIfADJB4ANAJgh8AMgEgQ8AmSDwASATBD5KxRuqgOqweBrGotcKl7yhCqgWgY/k+gU7b6gCqkVJB8n1e0sVb6gCqsUMH8n1e0sVb6gCqsUbrzAWvKUKqAZvvELpeEsVUD/U8AEgEwQ+AGSCwAeATBD4AJAJAh8AMkHgQxJr3AA5oC0TrHEDZGKkGb7tn7T9iu3vFX/3TAnb/2P7QvFnapRzIr1+SyEAaJdRSzpPS/pmROyX9M1iu5f/jIiPFX8eGfGcSIw1boA8jFrSOSLpvuLn5yR9W9LvjPidKBlr3AB5GDXwPxwR7xY//5ukD/c57sdsdyUtSno2Iv6m10G2j0k6Jkm7d+8ecWjYCJZCANpvzcC3/aqkO3rsOrFyIyLCdr+V2PZExBXb+yR9y/Y/R8S/rj4oIk5JOiUtL5625ugBAOu2ZuBHxAP99tn+vu07I+Jd23dKeq/Pd1wp/p6z/W1JPyfplsAHAIzPqDdtpyQ9Wfz8pKQXVh9ge6vt24uft0v6BUlvj3heAMAGjRr4z0p60Pb3JD1QbMt2x/bp4piPSura/o6k17RcwyfwAaBkI920jYhrku7v8XlX0tHi53+Q9LOjnAcAMDqWVgCATBD4AJAJAh9DY8E1oFlYPA1DYcE1oHmY4Zcs5ay4yhk2C64BzcMMv0QpZ8XDftfM/EKSNXNuLrh2Y3GJBdeAhiDwS9RrVjxs6A7zXSl/4bDgGtA8BH6JUs6Kh/mulL9wJBZcA5qGwC9RylnxMN9FGQbImyPquShlp9OJbrdb9TBaJ1UNH0A92Z6JiE6vfczwM0MZBsgXbZk1UeeHmOo8NgDrxwy/BlI/xJSybMMDVkB7EPg1kLJ7JnVAp+7sAVAdSjprKKOccbN7ZsIauXsm9ROwKccGoFrM8Acoq5yRsl0zdeslD1gB7UHgD1BmOSNV98w4AprOHqAdCPwBmvqgEgENoBcCfwDKGQDahMBfA7NlAG1Blw4AZILAB4BMEPgAkAkCv+ZYxwZAKty0rTHWsQGQEjP8GuNF4QBSIvBrjHVsAKRESadQxzdB8eAXgJQIfNW7Vs6DXwBSoaQjauUA8kDgi1o5gDy0sqSz0Xo8tXIAOWhd4A9bj09ZK6/jDWAAaF3gV/0O1jrfAAaQt9bV8Kuux3MDGEBdtW6GX2Y9vlfppqlvyQLQfo6I4f9j+1FJvyvpo5IORkS3z3GHJf2xpAlJpyPi2bW+u9PpRLfb8+tqYVDphho+gKrYnomITq99o87w35T0SUl/MeDkE5JOSnpQ0mVJ521PRcTbI567UoPuFfCwFIA6GqmGHxHvRMTFNQ47KGk2IuYi4rqkc5KOjHLeOqj6XgEAbFQZNfwdki6t2L4s6d5eB9o+JumYJO3evTv5QFKWWujdB9A0awa+7Vcl3dFj14mIeCHlYCLilKRT0nINP+V3j6NdktINgCZZM/Aj4oERz3FF0q4V2zuLz0pVdX8+AFStjD7885L2277L9hZJj0maKuG8H0DNHUDuRqrh2/6EpC9LmpT097YvRMTHbf+UltsvH46IRdvHJb2k5bbMMxHx1sgj3yBq7gByN1If/jjVvQ8fAOpoUB9+65ZWAAD0RuADQCYIfADIBIEPAJkg8AEgEwQ+AGSitm2Ztq9Kmi/hVNsl/aCE8zQd12n9uFbrw3Van41epz0RMdlrR20Dvyy2u/16VvH/uE7rx7VaH67T+qS8TpR0ACATBD4AZILAL5Zjxpq4TuvHtVofrtP6JLtO2dfwASAXzPABIBMEPgBkgsCXZPv3bL9h+4Ltl4v1/LGK7S/Z/m5xrf7a9oeqHlMd2X7U9lu2l2zTdriK7cO2L9qetf101eOpK9tnbL9n+81U30ngL/tSRNwTER+T9HeSPl/xeOrqFUk/ExH3SPoXSZ+reDx19aakT0p6veqB1I3tCUknJT0k6W5Jj9u+u9pR1dZXJR1O+YUEvqSI+I8Vmz8uiTvZPUTEyxGxWGxOa/n9xFglIt6JiItVj6OmDkqajYi5iLgu6ZykIxWPqZYi4nVJP0z5nSO94rBNbH9R0q9Jel/SL1Y8nCb4tKTnqx4EGmeHpEsrti9LureisWQnm8C3/aqkO3rsOhERL0TECUknbH9O0nFJz5Q6wJpY6zoVx5yQtCjpbJljq5P1XCegbrIJ/Ih4YJ2HnpX0ojIN/LWuk+2nJP2ypPsj44c4NvC/J3zQFUm7VmzvLD5DCajhS7K9f8XmEUnfrWosdWb7sKTflvRIRPyo6vGgkc5L2m/7LttbJD0maariMWWDJ20l2f4rSR+RtKTlJZk/ExHMOlaxPSvpdknXio+mI+IzFQ6plmx/QtKXJU1K+ndJFyLi45UOqkZsPyzpjyRNSDoTEV+sdkT1ZPvrku7T8vLI35f0TER8ZaTvJPABIA+UdAAgEwQ+AGSCwAeATBD4AJAJAh8AMkHgA0AmCHwAyMT/AgJb6ug6qPc6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 原始x散点图\n",
    "plt.plot(data['x'],data['y'],'.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 普通线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "\n",
    "def myplot(x, y, y_pred, sub, title):\n",
    "    plt.subplot(sub)\n",
    "    plt.tight_layout()\n",
    "    plt.plot(x, y_pred)\n",
    "    plt.plot(x, y, '.')\n",
    "    plt.title(title)\n",
    "\n",
    "\n",
    "def summary(y, y_pred, intercept_, coef_):\n",
    "    rss = sum((y_pred - y)**2)\n",
    "    ret = [rss]\n",
    "    ret.extend([intercept_])\n",
    "    ret.extend(coef_)\n",
    "    return ret\n",
    "\n",
    "\n",
    "def linear_regression(data, power, models_to_plot):\n",
    "    #设置预测变量x,x_2,x_3...\n",
    "    predictors = ['x']\n",
    "    if power >= 2:\n",
    "        predictors.extend(['x_%d' % i for i in range(2, power + 1)])\n",
    "\n",
    "    #线性拟合\n",
    "    linreg = LinearRegression(normalize=True)\n",
    "    linreg.fit(data[predictors], data['y'])\n",
    "    y_pred = linreg.predict(data[predictors])\n",
    "\n",
    "    # 绘制指定的幂的图形\n",
    "    if power in models_to_plot:\n",
    "        myplot(data['x'], data['y'], y_pred, \\\n",
    "               models_to_plot[power], 'power=%d' % power)\n",
    "\n",
    "    # 记录模型拟合效果rss、截距和系数\n",
    "    return summary(data['y'], y_pred, linreg.intercept_, linreg.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEYCAYAAACQgLsAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABE+0lEQVR4nO3dd3hUZfbA8e+ZmTRCaKEJJIRepBMCoqAoKmJBXQs2cF3F7vrT1bVsX7GsuvaGiw0LFkRAECkqBgFDkd5LIAGkJECAkDIz7++PGWQIk2SSTDJ3JufzPHlkMnMzJ/GeOW+77xVjDEoppVQw2UIdgFJKqcijxUUppVTQaXFRSikVdFpclFJKBZ0WF6WUUkGnxUUppVTQaXFRSikVdFpc1G9E5BoRWSAi+SLyQ6jjUSrciEgjEdknIvNDHUuoOUIdgKo+IiKAGGPcAR6SC7wIdAbOra64lLK6SuTOcc8A69CGu/4BqoOIZIrIoyKyVkQOiMi7IhLrfe42EdksIrkiMlVEWni//08RecX77ygROSoiz3ofx4lIgYg08j4e4O1hHBSRFSJyjs97/yAiY0XkJyAfaBto3MaYOcaYz4BdQfpTKFUh4Zo73uMHAt2Ad6v+lwh/Wlyqzw3AhUA7oCPwFxE5F3gKuAY4DdgOTPS+fh5wjvff/YBfgcHex2cAG4wxuSLSEpgOPAE0Av4ETBKRJj7vfRMwBkgAtovI695k8ve1sjp+eaWqIOxyR0TswKvAPYDuqYUWl+r0qjEmyxiTC4wFrsOTNO8YY5YZYwqBR4EzRCQFWAh0EJFEPIkxHmgpInWBs/EkEMCNwAxjzAxjjNsYMxtYAgz3ee/3jDFrjDFOY0yxMeYuY0yDUr561MDfQqmKCMfcuQ/42RiztNr+KmFGi0v1yfL593aghfdr+/FvGmOOADlAS2PMMTwn+tl4EmQesAA4k5MTpDVwtW8LCjgLT2vO33srFW7CKne8w3P3AY9X9NhIphP61SfJ59/JeOYxduE5wQEQkXggEdjp/dY8PBPpvYHF3scXAmnAj97XZAETjDG3lfHeJ3XLReRNPK02f7YbY04P4PdRqqaEW+6k4SlQaz3rAIgD4kTkVzzFz1XWLxuxjDH6FeQvIBNYBbTCM7Y7H3gSGArsA3oBMcBLwHyf4y4A8oC53senex+v8XlNEp4x5QsBOxCLZ7y5lff5H4BbKxn38Z93B56EjAWiQv331K/a8xWOueONp7nP1x+Bn4Hmof57hvJLh8Wqz8fALGArsAV4whgzB/grMAnYjWfCcqTPMQvwtHqOt7TWAgU+jzHGZAEjgMfwJFsW8BDBGeK8CTgGvAEM8v777SD8XKUqIqxyxxhTaIz59fgXcAgo9v671hJv5VVBJCKZeFpAc0Idi1LhRHMncmjPRSmlVNBpcVFKKRV0OiymlFIq6LTnopRSKuhCcp1L48aNTUpKSijeWim/li5dut8Y06T8V4aW5o6ykrLyJiTFJSUlhSVLloTirZXyS0S2l/+q0NPcUVZSVt7osJhSSqmg0+KiIp7T5WZPXkGow1C1RVYGpD/v+W8tpnuLqYhljOH7DXt5csZ64qPtTL7rTGw2CXVYKpJlZcD7l4GrCOzRMHoqJKWFOqqQ0OKiItKaXYd4csY6ftqcQ0piHe48pxOidUVVt8x0T2ExLs9/M9O1uCh1kqwMT2KkDAqr5Pj1UAHPzdrApGXZNIiL4h+XduX6/q2JdugIsKoBKYM8PZbjPZeUQf5fF6b5VRFaXNSpwrBrf7TQyVvztjAufStuN4wZ1Ja7hrSnflxUqENTtUlSmidfyiocYZhflaHFRZ2qIl37ELfAXG7DZ0uyeH7WRvYfKeSSHqfx52GdSWpUp8ZjUQrw5EFZuRDsoTOL9oK0uKhTVaRrH6IWmDGGeRv38dSM9WzYc5jU1g15e1Rfeic3rJH3DwYReQe4BNhrjOkW6nhUDSkvvypSLCzcC9Liok4VSNceQjZ5uXZXHk99s470TftpnViHN27ow7BuzZHwm7F/D3gV+CDEcaiaVFZ+BVosjhegQ9mWXUAQlOKiLbAIVF7XHspugVWmq17OMXvyCnh+1gY+X5pN/bgo/nZJV24cEL6T9caYH0UkJdRxqBAoLb8CabD5FiCbHWwOcFP2KEMIBKvn8h7aAqt9SmuBVaarXsYxRwudjPtxK+N+3IrLbfjDmW2499wO1K9TOybrRWQMMAYgOTk5xNGoCqlgI6s46Uwc9ihwgbFHkZOYRvSxYurFOk70zH0LkBvoOwrqJ0XmnIu2wGoxfy2wyiwI8NO9d7Xsx+dLsnh+9kb2HS7k4h6n8ecLO5OcWLsm640x44BxAKmpqXqPjHBRRoPpUH4xS3fksm73YdbuzmPzniPsPVzAgfxi+sgjDLCtY5G7C8s+yANmEe2w0aRuDMmN6nB+vVaMkijsAPZopOf1pw6tWWCCX+dcVPBVZkFAie79L7ZuPPJSOhv2HKZPcgPevLEvfVuHz2S9Ur6NLOMqYveK2Xy8tj7pm/ezKvsgbm8zIalRHJ2aJdCvTUOa1I0lsW43ou02UoBrMRwucLLvcCF7Dxeydf9Rnlldj69dngK0uU4vEhfHce7hPZzdsQnRu5dYZoK/xoqLdu3DWEVbQpVZEODt3u+1NeXN7S1Y/vVarqy7me7DLuaMsweG42S9qu1SBuG2R4ETioyde36qwwrZQq+kBtx7bgcGtkuka4t6JMRWbHjX6XKTmXMWGdsOsHnDXqYu38knGTtoFB/Ns83mcK6rCLHABH+NFRft2oepyi51rOCCAGOP4o3cNJ5bV5+BMVv5LO4p7E4nMn8StLPO8spgEpFPgHOAxiKSDfzdGDM+tFGpqjLGsHBLDu/8JOTme3oYR5sPYGT/83m3W/MqX9jrsNto3zSB9k0TuL5/MkVONz9t3s/nS7N4c91pDLTbiRYDtihscY2Q9OdDMkSmw2KqbNW53DgpjYLrJ7N03jRe3dacJRsacPPAFB6KX4Mj3WnJ5ZXBZIy5LtQxqOAxxjBr7R5enLOJdbvzSIyP5oYhF3N92p00rx/raagt+cz/8uMqzJFEO2wM6dyUIZ2bknu0O9/OTWbn8jlkHYvln18/TDROxBFT40NkwVqKrC2wSBXo/EkFudyGScuyeX7WUfbkncnw7s2ZfWFnUhrHQ9YRWPDfE+8Zl+jZwtxiq2GUghMX9P539kZWZh+iTeN4/vO7HlzWqwWxUXbPi0obAajIyEAARahRfDSXX3YlBReNYN1n/8C+yYngxuUs5PDa72gQbsVFW2ARLND5kwqYv2k/Y2esY93uPHolNeC16/uQmtLI/3vGJcLMRywxQalqodI+0L3f39mgL48trsPhTT8xPH4zfzxvGGefOxiHvcS1V6WNAAQ6MlDB4enYKDu9B1+K2fY2bmcRRcbOs+n7GLbrUQYMGUFUyoAg/YFKp8Niqny+8ydltZ7KaVlt+PUwT32zjh827KNVwzheua43l/Q4zf9k/fH3TH/eslcgqwhXRm/DvH8ZxllII+OgDTfzl7j3sbucyM+ToLOfD/7SRgACHRmozPB0UhoyeiqSmU6xJPD3uX/BllmM873/kTXiU9r2Pjewv0ElG5ZaXGqDQAsCVH431zKe23u4gBdmb+TTxVmeyfquv9Jr8CVEp7QoP/aSyadDZKqmlPKBvn3pt7R0FuLATbQ4ebT1Rhw7ypkjLG0EINCRgcoOT3sbafXSnwecIG7AyZeTJhJ7IIm7zmlf+g30qrhvmRaXSBdoQbDZAQG3s/QTqazWk5/njjXry9vpW3lz3haKnG4e73GYW7Y8gWwrhh3jAztZdYhMhUqJD/T8FgMZO3kV6zLq8nGMA7u4sNujsXe/Anb+XP4Hf2krKMtaWenb+KvK8LTP72K3RyHJg3lu1kYyMg/w4rW9aBQffeoxVVzMo8Ul0gVcENzeA0zpJ1JZraeTlhVH811BRx577nv25BUy7PTm/PmizrRZ9yZsLK74yVrVITKLXLGswoxPw2ZTXC9e+Hw5KUd+4bJe5yH9piI7F5w4p5p1Df455q9hOOjBKv8ukjKIB1r147SMLP4xdQ2XvJzOazf0OXVH8Sou5tHiEukCLAin9Fz8nUhldeG9z21f+i0vbWnGl3Nt9EyK49Xr+9DPvhnWvenpeVRl5VllTnYLb0murM/dsh9vb01kzuRpTIgaS4zDiWyaCgNLfNAHcl1XRQX7MgCfGAW4vn8y3VvW586PlnLNWwt59qqeXN675cmvr0JvSYtLpPN3gpTW1YbyT6RSkmjTnsM8NVf4bn0fWjaI4+XrOnNJ99Ow7Vx88of7sKfhWE7lWniVOdn1nuaqkg7mF/HHicuZt3EfL7bIIuaACzHumjuPqukyAF/dW9Xn63vP4o4Pl3L/p8vZkZvPvee2P7HIpgpFU4tLbVBytVdZXe0Knkj7DhfywpyNTMzYQXyMg0cv6szogSkn1veX/HA/llP5rn3J3yUQNZCgKgKUGDpdvfMQd3y4lL15hTxxeTdGtGiMfPBRzZ5H1XAZgD8NcpbzYaefeC2qOf+dvZEdufk8dWV3okoup64gLS61TZBa8seKXIyfv5U3fthCodPNqDNSuO+8DqdODIb6w72GElSFsRINrh8G/I/bf7DTsE40n94+wDsX0To051F1DLf58v7uDlcR99mj6dr1fn5ZPoVnc8/kgVtuPNFIrAQtLpEkkInrKn7Yu92Gyb/s5LlZG9h9qIALujbjkYs607ZJXf8H1OSHe2m/f3UnqApvPg0ul7OIjO+n0rv173n1+j40rhtz4nWRdB75udWFOAs5P/NZzotyU7hrMs+87ebhW0cRF125AqPFJVIEOnFdhQ/7BVv2M3b6OtbsyqNHq/q8eG0v+rdNLP/AmkjKkr9/VeZ2VO2SMghjj8btLKTI2KnXZQgTrutf5WEhyyrtVhciYNzYcBMjLursWsjodzoz/ubUCu/cDFpcIkdFhrsq+GG/ee8Rnv5mHXPW7aVlgzhevLYXl/VsUfrFVzXJ383GnIUw40EwRleIqXLtrd+D/9T5F81yl9DljOHcPnxEZN/iobQ7WfpcR2azR3PG4BG8NecA7/2Uyb3ndajw22hxiRTVMLex/0ghL87ZyCcZWdSJsvPnYZ35/ZkpVRqHDapyWmDU5MoeFZa27DvC6HcyyD2azGs3Xs6Qzk0j/7qokp8Vvney9Lle56ykND5vf4DuLetX6m20uESKIM5tFBS7GD9/G2/8sIVjxS5u7J/Mfed1INF3/NkKAmiB6QoxVZql2w/wh/cX47AJn445g+6t6teO66LKu17N5/EpF1ZWgBaXSBLoBpOlcLsNU1bs5NmZG9h1qIChXZrx6PDOtCttsj7UAmyBRdyHg6qypfO/Zd6sLxlYpxd/vm0UrRPjPU/UluuiamAeVItLJKpE62vR1hzGTl/Hqp2H6N6yPs9f04sz2gUwWR9KFWiBqQhXgcbU/O+m03fezfSyObE5v0Ly+0Ki95hQL52PIFpcrKQyY73+jqlA62vLviM8NWM9c9btoUX9WF64ticjera0xmR9ILSIqPI2Z/XJj08ydpD13RTOcDix4wZX8cn5oddFBY0WF6uozFhvaccE0PrKOVLIS3M38dHPO4iLsvPQhZ34w1ltrDNZr1SgSmtMlciPqb3e4NH50dyaMhDb/qml54c2WIJCi4tVlNXbKK1HU9oxZbS+CopdvLcgk9e+20x+sYvr0pK4f2jHky8WUyqclNaY8skPt7OI9Qu/4eLud/DwtRchu3VOrrppcbGK0hKkrB5NWT2UEq0vt9swbeUu/jNzAzsPHmNol6Y8clFn2jdNqMFfUqlqUFpjyntxpHEWUWjsRLUfzEsje3luQay9k2qnxcUqSkuQsno0AY4P/7w1hydnrGNF9iFOb1GPZ6/qwcD2jWvoF1OqBvgpFqZVPz7p8grZy2YT1/Fs7rvpeuzhMpcYAbS4WIm/1lR58ydltMC27jvC09+sZ9baPTSvF8vzV/fkit5hNFmvVGVkZWC2pTNxXzKPLa7DyH738OQV3fW8r2FaXKyuEqtXco8W8fLcTXy4aDsxDht/uqAjfzirbaU3oFMqbHiHkY2zkMuNg9yuL3HnFcM9hSXSr7y3GC0uVlUyEXyToZQkKSh28f6CTF79fjNHC52MTEvm/4Z2pEmCTtarWiIzHbezEBtuosXJnSm7TxSWSL/y3mK0uFhReev2SzxnWvVj2srd/GfmerIPHGNIpyY8OrwLHZvpZL2qXSbltmG4cRAtTmyOGKTNqSvHIvrKewvR4mJFZSVCieeyf5nFPVOLWZ51kC6n1ePDP/TgrA46Wa9qn/Hzt/HvhTHs6PQCf2y3x1NYAllZqaqFFhcrKisRvM8ZVxHFOLhvYTy7E47x7FU9uLJPK10No2qlCYu28++v1zK8e3PuHXkRtpL3YtEr72ucFhcrKiMRDjTqxZftXuTA2u9YKt04Z+hF3DZIJ+tV7fX5kiz++tVqhnZpxksje3uuY/FHr22pUVpcrKpEIhQ6XXywYDuvfLeJI4UJXNvvXl46vwNNE2JDGKRSofX1yl38edJKBnVozGs39I7cu0eGIS0uFmeMYfqq3Twzcz1Zucc4u2MTHhvehU7NdbI+EojIMOAlwA78zxjzdIhDChtz1u7h/onLSW3diHE3pRLj0N67lQSluGiCVI+l23N5Yvo6ftlxkM7NE/jgljQGd2wS6rBUkIiIHXgNOB/IBhaLyFRjzNrQRmZ9P23ez10fL6Nri3qMvzlVh4UtqMrFRRMk+LbnHOWZmeuZsepXmibE8J+revA7nayPRGnAZmPMVgARmQiMADR3yrB0+wFu+2AJbRLjef/3aSTERoU6JOVHMHoumiBBcjC/iJfnbmbCokyi7Db+b2hHbhvchjrROnoZoVoCWT6Ps4H+JV8kImOAMQDJyck1E5lVlLhgeM2uQ9z8bgZNE2KYcGsaDeOjQx2hKkUwPrU0QQJRxtYThU4XExZu55XvNnO4oJhrUpN44PyONK2nk/UKjDHjgHEAqampJsTh1JwSFwxnXzaRUVOKSYhx8OGt/XUxi8XVWJO41iYIlHrFvTGGb1b/ytPfrGdHbj6DOzbhseGd6dy8XqgjVjVjJ5Dk87iV93sKTrpg2LiKmDrlM0Su4MNb+9OqYZ1QR6fKEYzioglSHj9X3C91d2Ds9LUs807Wv39LGmfrZH1tsxjoICJt8OTMSOD60IZkIT4XDBcaO3td8cxIW0zTgvpAmm5EaXHBKC6aIOXxueLebY/ivxub8ur0BTRJiOF/57o4N3YxtrhYQItLbWKMcYrIPcC3eFZavmOMWRPisKwjKY28aybxxaRPyDwWyz+iPsC2uBiWvQzDnoaZj+hGlBZW5eKiCRKApDSOXPslGT9M4Y3tp+HIPMCETmtJ69KemDmPexLkR02Q2sgYMwOYEeo4rCivoJgbvjVsPHoxs1KXYFtRfKL3v26KbkRpcUGZc9EEKV2R082ERdt5ee4R8goG82CXg9y94wFkRzFkCRi350sTRKnf5Bc5+cN7i1m3O49xo/rSOj4RVr92oqfSZQRsX6gbUVqYrnGtJsYYZq7+ladnrmd7Tj5ntW/MY8O70HXL27DN2wIzNrDZANEEUcqr0Oni9glLWbr9AK9c14dzOzcDmp26316zrjrnYmFaXKrBLzsOMHb6OpZsP0DHZnV57/f9OLtjE0QEXCV2PB72NBzL0QRRCih2ubnn419I37SfZ6/qwcU9TjvxZMmNJ3UjSkvT4hJEWbn5PDNzPV+v3E3jujE8dWV3ru7b6uRdWnXrb6X8crrc3P/pcmav3cM/Lzudq1OTyj9IWZYWlyA4dKyY17/fzLs/ZWKzwb3ntuf2s9tRN6aUP6+2uJQ6idttePiLlUxfuZvHh3dh9MCUUIekqkiLSxUUOd18/PN2Xpq7iYPHirmydyv+dGFHTqsfF+rQlAobbrfh8a9W8eUvO3nw/I7c1mY/pE/Wnn2Y0+JSCcYYvl2zh2dmrmfb/qOc2T6Rx4Z34XTXBlj5+omJeR36UqpMxhj+OmU1GxbP5YOOexlcfx+8r9evRAItLhW0IusgY6evIyMzlw5N6/Luzf04p1MTJHvxiS1ebHZAwO3UBFGqFMYY/jZlDesy5vBp3FM4spyQrcvzI4UWlwBlH8jn2W83MGX5LhrXjWbsFd24NjXpxGT9SVu8uL1HGU0QpfwwxvDPaWuZsGg777Xfg2OnE9Hl+RFFi0s58gqKef37Lbzz0zYEuGdIe+44x89kvc8WL6f0XDRBlPqN223429TVfLhoB7ee1YazezRGPnhPl+dHGC0upSh2ufkkYwcvztnEgfwirujdkj9d0IkWDUqZrC+5xBh0zkWpEptLut2GxyavYuLiLG4/uy2PDOvsuf5Ll+dHHC0uJRhjmL12D09/s56t+49yRttEHr+4C93cG2DV62Wf/P4u8lKqtipxqwnXTVN46OcYvly2k3vPbc8D53f0FBbQ5fkRSIuLj5XZnsn6n7fl0rZJPP8blcp5XZqePFmvE/RKBabE/Vi+nvoZX+48lwfO78h953UIdXSqmmlxAXYePMazM9fz1fJdJMZH8+/LuzGyXxJRfifrdYJeqYD43I+lCAfv72rFXy7uwq2D2oY6MlUDanVxOVxQzOs/bGH8fM9k/V3ntOPOc9qREBt18gtTSuwHphP0SpXPe6uJr76ayOQDbRj5u6u4Rrd0qTVqZXFx+kzW5xz1TtZf2ImWgU7Wa69FqXJl5ebz+2nF7Dh4ES9f35th3ZqHOiRVg2pVcTHGMHfdXp76Zh1b9h2lf5tGvHdxV7q3ql/+waVNOOqtVpU6xarsQ9zy/mIKi1188Ic0BrRNPPVFmjsRLbKLi8/Ju9rWibHT17Fwaw5tm8Tz9qhUhh6frE+v5AleYjWMTvSriBDoh34pr1ucPpP5syfTO7YnD915Ex2aJfg/VnMnokVucfGevMZVRDEO/lbwKJl1uvGvEadzXVqyZ7K+rBM8kATTiX4VaQL90PfzOtOqH19OnczwZbfTx+bE5v4KKUoF/ByvuRPxIra4FG6eR5SzEBtubMZwX7s99LnxHur5TtaXdoIHmmA60a8iTaAf+iVeV7zlR/4v3UHymplERzmx4wZXcenHa+5EvIgrLk6Xm4mLs5g7P47XjYNocWKLiuacXp1h8csn90RKnuBxiZD+PBzKLjvBfHs1OtGvIklpOVHy/PZ5ndsexWs/55J8+E3STm+PbWuM/6JRcjRAcyeihWdx8TNkZYzh+w17eXLGejbvPUJaSl+yUyfS4dhyT4LM9LONt+8J7vsamx1sDnDjP0FK9moGPRiSP4NSFVbecG9pOVGyB5+Uhhk1hdULpvP5mnwec75NtMPpKSz+9gYrbTRAi0rECr/i4uckXWP3TNYv2JJDm8bxvHVTXy7o2sy7tcR5npZXaT2R4ye472vcQN9RUD/p1CTUsWIVrgId7vWXEyXO9UP5xfx9QTRfLe/PU01mEXPEhRzfJv9YzqkNLs2bWif8ikuJLSVmTPuce7LOoUFcFP+4tCs3DGh94sr64wIZ3y35mp7X61ixiiwV/YAv5Vz/ZtVu/jplDQfyi/i/oR25pkNjZMKnFcsvzZuIF37FJWUQxh6N21lEkdvO+7taMWZQW+4a0p76cVH+jwlkfDfQMWAdK1bhqqIf8CXO9az4bjwxYQnfrtnD6S3q8d7v+9GtpfcasWDll4oYYoyp8TdNTU01S5YsqfBxTpebz5ZkM3vWNLoUrMDRbjBXX34lSY3qVOwH6cVbqgQRWWqMSQ11HOWpbO78phLn/qH8Yl79fhPvL9hOb9nI/3XYS+o5l+JoPaDycaiIUFbehEXPxRjDDxv38dSMdWzcc4TU1r05/+Lr6Z3csOI/TC/eUhHov7M38uuhY9wzpAPJiWU0tiowib73cAEfLtrB+wsyySso5sHOB7kr60lsmcUwYbzmjiqT5YvL2l15PDljHflbFnBV3c30uOhi+g8+48R9ICpKJxZVBHK53Xy1fBeTlu3kyt4tuXtIe1Iax1e4p+JyG5Zk5vLpkiymrdiF0204r3Mz/nRhRzpvehsyizV3VECsVVx8EmFP/R489+0GvliWzaDYrXwW9xR2pxNJnwRtq9Bi0olFFYEe6prHHfaVfLavNf9ZIXy+NJvfNdnJ00f/gsN4brctfnoaTpebzJyjrN19mPmb9jFn3V5yjxZRJ9rODf1bM3pgCm0ax3teXKy5owJnneLis12LUxzcX/wXlro7cOtZbXgwbi2OdGf5LaZAWmk6sagijTd3ElxF/MEezRU3fs5ne1pQd/EsxFWMiBtncSFvvvsesxOLsQkUFrspcLrIPnCMIqebPrKRwdEbuCnlLDr2PY+zOzWhbkyJjwfNHVUBVSouInI18A+gC5BmjKn8TGNmOsZZiOBG3Ibrmm3nmetu9YwfZx2BBc+XfdVvWRd8laQXb6kQC3bu+A71NtqXwR1nPwhtb8C8PxG3sxhjdyApZ5FQ5EAEEuNtdC5eT9/6a6jfqBm91jyNuIuRnV/B0KkQc5r/99LcUQGqas9lNXAl8FaVI0kZhHF4lhjbHNFcNuJaOD4xWVqLyXdyXgSM2/Ol48HK+oKaO36Hq5LSkNHTkMx0bCmDuDspjbuPH5OVAe8/5Dlml+aOCr4qFRdjzDqg8pPrvpLSsI2ednIBKTnMVfKE922xGRvYbIDoeLCyvGDnzkmNLzh5P7DyNp7U3FHVoMbmXERkDDAGIDk52f+LfBMhkCXDJVts/vY0UirMVSh3Krujt+aOCrJyi4uIzAH83Z/0cWPMlEDfyBgzDhgHngvByj0gkCXDOsGoLCwkuRPoUnvNHVXNyi0uxpihNRHIKQJdMqwTjMqiQpI7FVlqr7mjqpF1liKXpC0rpSpO80ZZRFWXIl8BvAI0AaaLyHJjzIVBiQy0ZaUiVrXmjuaNsoCQbFwpIvuA7RU8rDGwvxrCqSyrxQPWi8lq8UDpMbU2xjSp6WAqqhK5E07/D0LFavGA9WKqcN6EpLhUhogssdKutVaLB6wXk9XiAWvGVJ2s+PtaLSarxQPWi6ky8djKf4lSSilVMVpclFJKBV04FZdxoQ6gBKvFA9aLyWrxgDVjqk5W/H2tFpPV4gHrxVTheMJmzkUppVT4CKeei1JKqTChxUUppVTQhU1xEZF/i8hKEVkuIrNEpIUFYnpWRNZ745osIg1CHM/VIrJGRNwiEtJljCIyTEQ2iMhmEXkkxLG8IyJ7RWR1KOMIFavljtXyxhuTJXLHSnnjjafSuRM2xQV41hjTwxjTC/ga+FuI4wGYDXQzxvQANgKPhjie4/cI+TGUQYiIHXgNuAjoClwnIl1DGNJ7wLAQvn+oWS13rJY3YIHcsWDeQBVyJ2yKizEmz+dhPBDylQjGmFnGGKf34SKgVYjjWWeM2RDKGLzSgM3GmK3GmCJgIjAiVMEYY34EckP1/qFmtdyxWt6AZXLHUnkDVcsd625c6YeIjAVGAYeAISEOp6RbgE9DHYRFtASyfB5nA/1DFIvC0rmjeXNCROWNpYpLefe/MMY8DjwuIo8C9wB/D3VM3tc8DjiBj6wQj6p9rJY7VsubQGNSwWOp4lKB+198BMygBopLeTGJyM3AJcB5pgYuGgrZ/XUqZieQ5PO4lfd7qppYLXesljeBxGQBEZU3YTPnIiIdfB6OANaHKpbjRGQY8DBwmTEmP9TxWMhioIOItBGRaGAkMDXEMdVaVssdzZtSRVTehM0V+iIyCegEuPFsOX6HMSakVV1ENgMxQI73W4uMMXeEMB7fe4QcBIJ7f52KxTIceBGwA+8YY8aGIg5vLJ8A5+DZNnwP8HdjzPhQxVPTrJY7VssbsE7uWClvvPFUOnfCprgopZQKH2EzLKaUUip8aHFRSikVdFpclFJKBZ0WF6WUUkGnxUUppVTQaXFRSikVdFpclFJKBZ0WF6WUUkGnxUUppVTQaXFRSikVdFpclFJKBZ0WF6WUUkGnxUUppVTQaXFRvxGRRiLyqYjkiMh+EflIROqFOi6lrEZErhGRBSKSLyI/lHiuo4hMEZF9IpIrIt+KSKcQhRoyWlwimHhU5P/xE0BDoA3QDmgG/KMaQlPKUiqRK7l47rvytJ/nGuC5yVcnPDmUAdS62yhrcakGIpIpIo+KyFoROSAi74pIrPe520Rks7dFM1VEWni//08RecX77ygROSoiz3ofx4lIgYg08j4e4G01HRSRFSJyjs97/yAiY0XkJyAfaFuB0NsAXxlj8owxh4DJwOlV/4so5V+45ooxZo4x5jNgl5/nMowx440xucaYYuAFoJOIJFbyzxSWtLhUnxuAC/H0ADoCfxGRc4GngGuA0/DcFXCi9/Xz8NzxDaAf8Csw2Pv4DGCDMSZXRFoC0/H0MhoBfwImiUgTn/e+CRgDJADbReR1b3L5+1rpc9xrwCUi0lBEGgK/A74J1h9EqVKEY65UxGDgV2NMTrmvjCBaXKrPq8aYLGNMLjAWuA5PEr1jjFlmjCkEHgXOEJEUYCGe+2cn4jkZxwMtRaQucDaehAK4EZhhjJlhjHEbY2YDS4DhPu/9njFmjTHGaYwpNsbcZYxpUMpXD5/jlgHReG4/mwO4gNer5a+j1AnhmCsBEZFWeBptD1Ti7xLWtLhUnyyff28HWni/th//pjHmCJ4P8ZbGmGN4Tvyz8STMPGABcCYnJ0xr4GrfFhVwFp7Wnb/3rojPgI14WnH1gC3Ah5X8WUoFKhxzpVzeHtIs4HVjzCfV9T5W5Qh1ABEsyeffyXjGZnfhOeEBEJF4IBHY6f3WPOBcoDew2Pv4QiAN+NH7mixggjHmtjLe2/g+EJE38bTi/NlujDk+r9ILuNsYc9TnuPllvI9SwRCOuVIm77DyLGCqMWZsIMdEHGOMfgX5C8gEVgGt8Iz1zgeeBIYC+/B8iMcALwHzfY67AMgD5nofn+59vMbnNUl4xpgvBOxALJ7x51be538Abq1k3N8DrwBx3q/XgQWh/nvqV+R+hXGuHP95d+ApZrFAlPe5enhWiL0a6r9vKL90WKz6fIyn5bIVz/DSE8aYOcBfgUnAbjwTmCN9jlmA50P9eMtrLVDg8xhjTBYwAngMT/JlAQ8RnCHOW4AUIBtPC7EtMDoIP1epsoRjrtwEHAPeAAZ5//2297kr8Cw0+L2IHPH5Sg7C+4YN8VZaFUQikomnRTQn1LEoZWWaK5FLey5KKaWCTouLUkqpoNNhMaWUUkGnPRellFJBF5LrXBo3bmxSUlJC8dZK+bV06dL9xpgm5b8ytDR3lJWUlTchKS4pKSksWbIkFG+tlF8isr38V4We5o6ykrLyRofFlH9ZGZD+vOe/SqnAae4Auv2L8icrA96/DFxFYI+G0VMhKS3UUSllfZo7v9GeizpVZronOYzL89/M9FBHpFR40Nz5jRYXdaqUQZ5Wl9g9/00ZVPprdQhAqRMCzZ1akDc6LKZOlZTm6c5npnuSo7RuvQ4BKHWyQHKnluSNFhflX1Ja+Se8vyGACEwSpSqkvNypJXkTlGExEXlHRPaKyOpg/DwVJsoaAqhMt78WDBX40ryppcobOqtoHlg0b4LVc3kPeBX4IEg/T4WD0oYAKtPtryVDBSW8h+ZN7VPW0FlF88DCeROU4mKM+dF7b2tV2/gbAqhItz8rw/P8oexaMVTgS/OmFitt6CzQ3AmDvNE5FxV8x7v9x1tTZa2YOd7qstnB5gA35a9QUypSBZI7YZI3NVZcRGQMMAYgOblW3ZAt/B1vJZWxcmzv4QLW7Mxj2/6j7DxYl+hmz5Cct4yVjh6smVKETX6iSUIMzerF0KJBHKmtG9Fnx484jre63EDfUVA/6USCpD9f9mq1WkJzJ4z55I5p1Q8AESn99YGsNvPt3ZTMm6S0gPK1JtRYcTHGjAPGAaSmpuo+/+HCz5iuadWPTXuP8NPm/SzcksPK7EP8mlfw2yFxUXZaNUyhYd2OxDhsJNoEp9uwIyefxZm5HMwvBqC/I4oPohxEAWKPRnpefyI5LDqOHAqaO+Fl3+FCFm7NYdeqedy8+T4cOCk2Dm4ofoy19s60aliHpIZxtE6M5/QW9ejWsj7tm9Ylyu5dX1XeajOf3o2xR5Pf5RryGvemSd0YHBbKHR0WU2XzaSW5XUXMmv4Ff8vJY+/hQgCSG9VhQNtGdG/VgO7eJGlYJ6rM1tnB/CIytuWyaGsbHlwXQ/KhZayW7nRelcAt9QpoXkuWaqrIsm3/UV7/fjOTf9mJ0224P2YeUeLEjhsRF/e338u8pheSfSCfrNxj/Lwtl/wiFwDRDhvtmtSlU7O6dGiWQJOEGBrERVE/LgqX23C40MnhAid78grYkRNLVP2xJOUt4/uCjix6+wDwHXab8HD8dG51FmLHjXEVISHMnaAUFxH5BDgHaCwi2cDfjTHjg/GzVei43YbVUd3pjAObMRQbOx/sakW/To0Y3LExA9s1JqlRnQr/3AZ1orng9OZccHpzzCVdWLQ1l8xFmYyfv40PF23nX33a8Dt7NHK89RWXGJFDZJo3kWHv4QKenL6OqSt2EWW3ceOA1lzZpyVdXYnYJ3wFriLs9mgGn38Fg5O6/nacy23Ytv8oa3YdYvXOQ2zcc4SMbbl8tXxXme/XuG4MrRO7cLhdH3rUi+W8ujHUibGz+2ABR3aeQfGOzzDGiQsHewriSA5R7oTkTpSpqalGtw23rj15BXz88w6+WJrNzoPHODNmKzc020HTHufRrf/5xEbZq+V9d+Tk88T0tcxau4fhDXbwUKd9tElKhpmPVHs3X0SWGmNSg/6Dg0xzx1rmrtvDw1+s5HChk5sHpnDroDY0TYg98YJKzH8cLXSSe7SIQ8eKOXSsGJsICbEOEmIdJNaNoW5M2X0C5/ZFrJo/nelbCnnQ9S4x4sTmiKmW3Ckrb3RYTP1m6fYDvPvTNmau/hWXMQzq0IQ/X9SZ87sMIy7aW1DKSpYqTiQmJ9Zh3KhUfty4j39Mi2fIz8m8l5PO2a4iRIfIlIUUFLt4YvpaPly0g87NE/hkzAA6Nks49YW+8yel5UeJ78fHOIiPcZBUydgcrQfQu/UAun7/LI55Tmy4cTuLsNVw7mhxqeXMjp/ZsWwW72S3ZNXOQwyO2cDjPYYwZOglpDSOP/nFZU0WBjqRGEABGtyxCTPuG8S/v17LyxnNOCPGTrR4Jv0jdYhMhY8DR4v4z/8m0HDvz/yj93lc97thxOxeCumV2E+sIhPwFWy8xbQ/G/PT87ichRQZO6v3Cv1qMHe0uNQGpZyUG5bMIeXr62lpinlMbDhibdhwIZu/gjPbABW4ODKQSfgKJFJslJ2xV3RnervG/GGSg36s5eK+XWlfA0NkSv2mRO5k5ebzzP8+4NmjfyUmyolt01RYfqz8odvS8qMiF01WdBVYUhoyeipm6498vDyP61c+hbsah8hK0i33I93xk/K7sZ7/ZmWQlZvP3R8vY8rkz7BTjEPcROPCbopPHn4qqaw9kQLZarwS97q4uMdpPHnfH5hefySTf1qF21lY8XtlWHTvJWVxJXJn6y/fccXrC+iQv4JY8Qw34SqCdVPKP69Ly49At+iv7H1iktJwnP0nbuqZQLQ3ZhPg8cYYxn00kZ3TnqhU7mjPJdL5nJTGVcTCuV9x86YcbDb4e58LsK+bCq4ixGYHBNzO0k/ysi7wKm+/pMx0z5BWIFful5CcWIdJdw7klfczKdw12ZMk9mgkkOMttO5fhZkSuTNtyudEx1zFlZdei0ybfOKc6jICti8s+7wuLT/Ku2iyirlzXHS7s3HPfw5ncREucRDd+izKuJQTgPWL53LTxnuJESeseLXCuaPFJdJ5W0bGVUShsfPchiaMabef21vvIqHzEEjzObGh/DHdsi7w8vdcyQ/3YU/DsZwKj/smxEbx59tG8fGkOuxePpsDDfvzf/V70LS8AzPTPev9dUGAqiif3Ckwdo456jGjz2IaNDr31ILQrGvlc6e07wcpd46/h230NFamT+Nfqxpxxa7m3FTOZg/bl31LB3x6aBXMHS0uEc7ZIpVPO7/CzuWz2Rjbi0fObUvajzfDziL4+QVPkgx68MQBwf7gLdmdP5Zz8vtVgN0m3HT11UzreCYPf7GSua/M540b+9K3dcNSjzna4gwcxoEdgz3Q3o5SAElpbLnoY76e9hkFjvo8zHvYFhXD4hf9542FcweApDR6XteP+GMZPDl9HYPaNz510Y5XQbGLD3cnc64tCihjNKMMOucSwXYePMbIcYt4fEkd9va8m//+6TbSZE3N3uO7IrdMDtClPVsw+e6BxEbZGTluIU9/s57DBcWnvC7nSCEjZ7i5sfgxNna9D9EhMVUBa3flceXXTj6Pu4Y70hpicxXXXN5AteSOiPCfq3rgsAsPfr4Cp8vt93Wz1u5hfmFb1l/wIZz7eKWGk7XnEqG+X7+X+z9djtPl5sVre3F575aeJwLdsThYAr1lcgV1bl6Pafecxb++Xsub87bwxdIsnkw9RiprMa3P5GDj3oz5YAnZB47xxk3X0bVzs6C8r6odNu05zI3jfyY+2s4ntw2g/tG6kPFizeUNVFvunFY/jicu78YfJy7npbmbeLDLoVPe48tl2bSoH0u3/ueC7YJKvY9eoR9JsjIw29KZlteWP/4UTZfm9Xjthj608Xe9igV2TQ2WldkH+fTLSfwl5xGicFKMg38W30TzqHwuGH4VXfsPLfdn6BX6tZxPTmyN7cq14xYhwKe3n3EifyIsbx7+YgWbl33H53FPY3cX/zavc+TgHkZ/F82Asy/ioQs7l/kz9Ar92iArA/P+ZbidhZxvHNzZ9nnuGX0hdaL9/C+ujvHhUMnKoEdmOt3bZEGuE8GNDSdPxryPYJDZX0ELHQ5TZfCZOHfbo3hK/o7b3YGJYwac3DCLsLx5ovE8ZsYvwRQXAW5wFsKMB6njdvNhlIOcFr2BsotLWbS4RIjCzfNweHdDjRYXf+q0D5u/whJJfD4UxGYHexS4ndhEwLg9X7pCTJXHd+dvp6G7rOKBMTfSwd92LpHAmzfRriIuETvFxoYIiNgQtxsbbqLESauDS/Hsq1o5Ef7pUzsczC/iyRUN+adxECMu7I5oaFMLVkWVdtOkuMSTr5jWFWKqLCmDcNujcDsNxTi4+NKraXdavVBHVX188sYGZKf8jq+22tjniufvUROIquTqsJK0uIS5vXkF3DQ+g205Lbnqggme1WC15U6OJRcnHL/ZGAR23YFSQHbdboy1/4NOzhVceuk1tGtSN7Jzp0TetBt6G/ec1pdNe46wYO0Q6uxcSJ/Bl+Ko4u+uE/phbE9eAde+tZC9hwv536hUBrZv7HmiNl2VHqRJVp3Qr50y9x/l+rcXcaTQyQd/6E8vNtaO3KmBvNGeS5jaf6SQ699eRIvDq/ioXx4tYxoB3uJSm+7kGEmTrKrmZGWwf/Vc/r2kHgV05JMxAzi9RX3Pzsa1IXdqIG+0uFhJgK2Jg/lF3Pi/n2lycAUTop/Evqz45L1/avpaFqVCrSIt8awM3O9dSgNXEa/hYO8Vn5Pcor7nOc2doNHiYhUBDmUdKXQy6p0Mtu4/yrd9D2FfWXxqK6uaLr5SypIqOAycueRbWjmLcIgbu7hIzvNZFaW5EzRaXKyirKEsb6vMmXwmd84R1uzK460b+9ImIRHWvO6/laXDRaq2CCB3jheKqSt2MWFpAhMcDuzi8tyArmTvRHMnKLS4WEVp3XFvq8y4inDj4GjBozx15VUM7doMaKatLKXKyR1cRRh7NNN6vcF986NJa9Mf53lTiN29UPOmGmlxsYrSuuPeVpkYFzZjeLDjPs7sl3TycZocqjYrJ3c8F0cWsn7hN1zU7XZeuLYXsVF2aD8wtHFHOC0uVuKvUKQMwilR4Da4bVEMPG9EaGJTyspKyR1jj8btLKTIOGjWYyh/uroPNlt5t8lSwaDFxeJ+Lm7Hc0WPcnViJldeORJJ7h/qkJQKC/sa9OQ/8f+mac5iup15MaMv0oZZTdLiYlVZGRxc+x1v/FyXnAa9uPCOu3HERf32nM6zKFWKrAxy1szlr7/UZ15+W165/jLPHKXmTY3S4mJFWRmY9y8lwVnEGzjIvfBz6vsWltpwBbFSlZGVgeu9S6nvKuJFHGSNmEiH44VF86ZG6Z0oLci9LR23swg7bmLERcuDS0886W/ZpVIKgC2LZ2KcRTi8udPh2HLPE5o3NU6LiwV9kZNCkXHgxo6t5Dr8arj1qVKRYPrK3TyyrB5OcWDEfvI1LJo3NU6HxSzmu/V7+HNGLLs6vsAf2+3xbJ3v233XK4iVOsVnS7J4ZNJK+iSn4bpgCrKrxDUsmjc1TouLhWzbf5Q/TlxOl+b1uP2GYUi03f8L9doWpX4zaWk2D3+xkkEdGvPWTX09d19t5+caFs2bGqXFxSLyi5zcMWEpdpvw1k19iSutsCilfvPNqt089MUKzmyfyNujUj0XRypLCMqci4gME5ENIrJZRB4Jxs+sTYwxPPrlKjbuPczLI3uT1KhOqENSNURzp/J+2LCX+yb+Qu/khoy7SQuL1VS5uIiIHXgNuAjoClwnIl2r+nNrkw8WbmfK8l08MLQjgzs2CXU4qoZo7lTeLzsOcMeHS+nYLIF3bu5HfIwOwlhNMHouacBmY8xWY0wRMBHQS2EDtHT7AZ6YvpbzOjfl7iHtQx2OqlmaO5Ww8+AxbvtgKU0TYnn/lrQT14ApSwlGcWkJZPk8zvZ+T/nKyvDclzsr47dv7T9SyN0fLeO0+nH895peuudR7aO5U54SeXO00Mmt7y+hsNjF+NGpNK4bE+IAVWlqrC8pImOAMQDJyck19bbW4OfqYGeLVO79+BcO5Bfx5V0DqV9HW1/Kv1qbOyXyxn3TFO6fZ2fDr3m8+/s0OjRLCHWEqgzB6LnsBHz2gKeV93snMcaMM8akGmNSmzSpZfMKfq4OfnbWBhZuzWHsFd099+5WtZHmTllK5M1P333F7LV7+NslXTlb5yYtLxjFZTHQQUTaiEg0MBKYGoSfGzlKXB28yN2Vt+Zt5Yb+yVzVt1Woo1Oho7lTFp+8cdmieGFTU67q24rRA1NCHZkKQJWHxYwxThG5B/gWsAPvGGPWVDmySOJzdXB2g7688sUq/pW4iZG9r9OdWmsxzZ1yePPm8PrvuXdhPMkNo3mqyWwku9DznOaOpQVlzsUYMwOYEYyfFbGS0jjcpDfPvDye8fIvYvJdyIRPAAG3U3dqraU0d8rmbJHKLdOdOJxL+G/hX7HNK4b5z8Gwp2HmI7rLsYXpxpU1xO02PPDZCpLzlhGDEzEucBXrTq1KleG/szeyOPMAf+2Wg81VfCJX1k3R3LE4LS415JXvNjN77R66nDEcccR451+idKdWpUqxcEsOb8zbwrWpSXQ94+KTc6XLCM0di9PLWmvAt2t+5YU5G7myT0suHt4TerQ4MVYMOm6sVAmH8ot54LPlpCTG87dLu0KM49RdjZt11dyxMC0u1Wz1zkPcP3E5PZMa8OQV3RGRU3dn1cRQ6jfGGB6bvIp9hwuZdOfAE1u7+MsbzR3L0mGxarQ3r4DbPlhCwzpRvD2qr26sp1QAvliazfRVu3nggo70TGoQ6nBUJWlxqSaF2xYy882HSTm2hv+N7kfThNhQh6SU5WXuP8rkqZN5qsksbm+TE+pwVBXosFgwedfdO2MbIjMe4Xp3MTdERWN3DcCzR6FSyq+sDFxbf2T64oOMt71F7BEXMuFTXWIcxrS4BIt3HyTjKgIDduPGLgbcxZ5JR00Qpfzz5o44CxljBIcYxLhPLDHW3AlLWlyCxbsPkhgXGAGbHU+V0WWSSpUpMx3jLMSGGyOC2OyeHNLcCWtaXIIlZRBOcYDb4LZFEXXxM3AsV5dJKlWOw80HEIUDB07sjhjP1ffHcjR3wpwWlyB5c2sis449yuiW2Vxy6TVI6/6hDkkpy3O7DffNj+KI83FeGnCUFj3P14ISIbS4BKqMTfJe+34zz367gUt7nsPF1/TEbtdFeEoB5W4u+ca8LXy/YR//HnE5Lc5Iqfn4VLXR4hIIPzf7IikNYwwvzNnEy3M3cUXvljx7VQ8cWliU8iglb477afN+np+1gRG9WnDjgNYhDFRVB/0kDISfm30VOl08+NkKXp67iWtSW/Hc1T21sCjly0/eHJeVm899n/xC2yZ1T+xcoSKK9lwCcfymRd4WWF6zAdz6vwwyMnP50wUduXtIe00OpUoqkTfHV37tO1zITeN/xuk2vHlj3xPbu6iIov9XA+Fzs6+Ncb249asifs0r4OXrenNZzxahjk4pa/LJm+NzLocLirn53Qz25BXy4a39ad+0bqijVNVEi0uAXC378camhrzwzSaa14tl4pgB9EluqHfDU6osPptLFhS7GPPBUjb8epi3R6fS17YJ0jV3IlVkF5dAPvgDeM2mPYd5fPJqMjJzubRnC564vBv146LKnbBUKiwF2mCqQMNq58Fj3DFhKat2HuLFa3sxpE6m5k6Ei9ziEsgHf1mvycqgYNM8Pvw1iadW1yM+2s5zV/fkd31anphf8TdhqQmiwlmgDaZycse36Py8NYe7PlpGodPN26NSOb9rM0ifqLkT4SK3uATywV/Kaw5v+onYT67A4S7mBuPA2fVlrr78ShLrxpx8fCkTlkqFrUAbTKW9zqfoGHsUEzq+wr+W1yW5UR3GjUo9MceiuRPxIre4lDx54xIh/fmTu/AlXpPdoC/jp62h3uIPuFeKcYgbu83FHa13gW9h8W2Zlbw7nlLhLJC88fc6b3Fwbf0Rm3ePPVex4deVc7i05138c8Tp1Nv3y8lzLJo7ES1yi4vvyRuXCDMfObULn5RG7lWfs3XJt3y2rzWbP/6FgfZ1tGjVCtu+aHAXIyVbVf6GAwY9GLrfU6lgKpE3xps3blsUX/V4gyWu9jhdBoihdfLztDmynE1xPXHO2kCTnA/ZmBfFX+x2ovDssXfDVdfTsnuv0ofRtKhErPAsLoFOJB4/edOf/60Lb1xF7F4xm0/XN+C79XtZtTMPOIPfNdnJZ3FPYTdOZH80XPSM/83zdJ5FhbNAcicpjYLmfdnw+T/p5izEjhvjNGzJmMnsuKuJttswQLppjt12ET3zNvDfgr8RhRN3dBSrezxG+/giEjoPoeXx99C8qXXCr7hUYoVWTuM06ksUYgxFxs49P9XhFzbRJ7khD13YiQu6NqPDxnHwnfPEyX8sx3+PRMeKVbgKIHfcbsN7CzJ55btNtDnWkI9jHAgubI4o7h51Mw+3G3jqz03/Bb5zgXFjN056J7ph0J9Pfo3mTa0TfsUlgBaQMYY1u/L4ZvVuvl+/j7W78+gjjzC0ziZcyWcyuuc5vN2+8ckT9EUBnvw6VqzCVTm5szevgAc/X0H6pv0M6tCYu4eMIsYxENk+H1IGEV/auR5I4dC8qXXCr7iUcSLvyMnni6VZfL1yN1v3H8VuE/q2bsifh3VmSOdBdGqW4FlGnJUBv0w8+SSvyMmvY8UqHJWRO9+v38uDn68gv8jJE5d344b+yd4l94mQ7HP7CH/DaoHmjuZNrRJ+xaXEiexskcp3a37lw5938OPGfdgEBrRN5NZBbRnWrTmN4qNPPr6soQE9+VUk81MEjDG8/sMWnpu1gc7N6/HKdb1o3zTB//GaO6oCwqe4lGgxFTTvy+dLs3n7k3nsyM2neb1Y7h/agZH9kmleP7b0n6MTi6o2KdnT8CkCx4pcPPTFCr5euZsRvVrwzO96EBtlL/1nae6oCgiP4lLiwqxJ3d7g6dUJ7D9SRK+kBjw2vDNDuzQLbMt7nVhUtUUZPY2l23N5fPJqNuw5zCMXdeb2wW3L39lbc0dVgLWKS2nLJDPTMT4XZm1dPJPT293Knee0o3+bRhXb7l4nFlUk8uaOM/lMnlxZj7nr9/DXet9wnrMQwf1bT2NP/R48/c16Jv+yk+b1YnlndD+GdG4a2Hto7qgKsE5xycrA+LSyim+czJEmfZmzdg/rVjXiYfeJC7OuumokD/cuZb8jnZBXtY23h2JcRbhwsLzgUZq27s/bWS050+4gSpw4jZ1bZkexcMZcomw27hnSnjvPaXfiXiqaOyrIqlRcRORq4B9AFyDNGLOk0j8sMx3jLMSGG2dxIS++/S6vu3IBaJ2YQvvurzM8YQsNu55LW9+T+3hSlHYVvlIWFPTc8fbs7cbwz54H6T5yIEcK08iY3578jT+wPaEPXRv0JC3WwRW9W9I6MV5zR1WrqvZcVgNXAm9VOZKUQeCIxu0sxtgdtO09jBdthfS3raN5j6FI8tWnHuM7piwCxu350slGZX1BzZ0iHNgx2BzRdG/fBtKfp27KIM4+72I47+JTj9HcUdWsSsXFGLMOCM4tfpPSsI2eBpnp2FIGcRWcOPmXv+K/NeW7esXYwGYDRCcbleUFO3eOXDOJqKwF1GvULLBeiOaOqmY1NuciImOAMQDJycn+X+Q7nuuzH1ipramSq1eGPe1/PzClwlgguZPYZRB0GRRY3oDmjqp25RYXEZkDNPfz1OPGmCmBvpExZhwwDiA1NdWUe4BuKaHCXEhyJ9Dlwpo7qpqVW1yMMUNrIpBT6JYSKsyFJHd0GyNlEdZZiuyPnvxKVZzmjbKAAC5pL52IXCEi2cAZwHQR+TY4YSkV2TR3VKQTY8qf/gj6m4rsA7ZX8LDGwP5qCKeyrBYPWC8mq8UDpcfU2hjTpKaDqahK5E44/T8IFavFA9aLqcJ5E5LiUhkissQYkxrqOI6zWjxgvZisFg9YM6bqZMXf12oxWS0esF5MlYmnSsNiSimllD9aXJRSSgVdOBWXcaEOoASrxQPWi8lq8YA1Y6pOVvx9rRaT1eIB68VU4XjCZs5FKaVU+AinnotSSqkwocVFKaVU0IVNcRGRf4vIShFZLiKzRKSFBWJ6VkTWe+OaLCINQhzP1SKyRkTcIhLSZYwiMkxENojIZhF5JMSxvCMie0VkdSjjCBWr5Y7V8sYbkyVyx0p5442n0rkTNsUFeNYY08MY0wv4GvhbiOMBmA10M8b0ADYCj4Y4nuP3CPkxlEGIiB14DbgI6ApcJyJdQxjSe8CwEL5/qFktd6yWN2CB3LFg3kAVcidsiosxJs/nYTwQ8pUIxphZxhin9+EioFWI41lnjNkQyhi80oDNxpitxpgiYCIwIlTBGGN+BHJD9f6hZrXcsVregGVyx1J5A1XLHWtvXFmCiIwFRgGHgCEhDqekW4BPQx2ERbQEsnweZwP9QxSLwtK5o3lzQkTljaWKS3n3vzDGPA48LiKPAvcAfw91TN7XPA44gY+sEI+qfayWO1bLm0BjUsFjqeJSgftffATMoAaKS3kxicjNwCXAeaYGLhoK2f11KmYnkOTzuJX3e6qaWC13rJY3gcRkARGVN2Ez5yIiHXwejgDWhyqW40RkGPAwcJkxJj/U8VjIYqCDiLQRkWhgJDA1xDHVWlbLHc2bUkVU3oTNFfoiMgnoBLjxbDl+hzEmpFVdRDYDMUCO91uLjDF3hDCeK4BXgCbAQWC5MebCEMUyHHgRsAPvGGPGhiIObyyfAOfg2TZ8D/B3Y8z4UMVT06yWO1bLG7BO7lgpb7zxVDp3wqa4KKWUCh9hMyymlFIqfGhxUUopFXRaXJRSSgWdFhellFJBp8VFKaVU0GlxUUopFXRaXJRSSgXd/wOXGJqsG6yR8QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "col = ['rss', 'intercept'] + ['coef_x_%d' % i for i in range(1, pow_max)]\n",
    "ind = ['pow_%d' % i for i in range(1, pow_max)]\n",
    "coef_matrix_linear = pd.DataFrame(index=ind, columns=col)\n",
    "\n",
    "# 设置显示4个图形，分别幂为，1，4，8，12\n",
    "models_to_plot = {1: 221, 4: 222, 8: 223, 12: 224}\n",
    "\n",
    "# 拟合所有幂的变量\n",
    "for i in range(1, pow_max):\n",
    "    coef_matrix_linear.iloc[i - 1, 0:i + 2] = linear_regression(\n",
    "        data, power=i, models_to_plot=models_to_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rss</th>\n",
       "      <th>intercept</th>\n",
       "      <th>coef_x_1</th>\n",
       "      <th>coef_x_2</th>\n",
       "      <th>coef_x_3</th>\n",
       "      <th>coef_x_4</th>\n",
       "      <th>coef_x_5</th>\n",
       "      <th>coef_x_6</th>\n",
       "      <th>coef_x_7</th>\n",
       "      <th>coef_x_8</th>\n",
       "      <th>coef_x_9</th>\n",
       "      <th>coef_x_10</th>\n",
       "      <th>coef_x_11</th>\n",
       "      <th>coef_x_12</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>pow_1</th>\n",
       "      <td>4.1</td>\n",
       "      <td>0.75</td>\n",
       "      <td>0.54</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_2</th>\n",
       "      <td>2.4</td>\n",
       "      <td>0.78</td>\n",
       "      <td>0.23</td>\n",
       "      <td>-0.14</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_3</th>\n",
       "      <td>0.89</td>\n",
       "      <td>0.99</td>\n",
       "      <td>0.11</td>\n",
       "      <td>-0.56</td>\n",
       "      <td>-0.13</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_4</th>\n",
       "      <td>0.85</td>\n",
       "      <td>0.96</td>\n",
       "      <td>0.043</td>\n",
       "      <td>-0.49</td>\n",
       "      <td>-0.038</td>\n",
       "      <td>0.021</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_5</th>\n",
       "      <td>0.75</td>\n",
       "      <td>0.96</td>\n",
       "      <td>-0.13</td>\n",
       "      <td>-0.58</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.18</td>\n",
       "      <td>0.03</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_6</th>\n",
       "      <td>0.74</td>\n",
       "      <td>0.95</td>\n",
       "      <td>-0.11</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.11</td>\n",
       "      <td>-0.0093</td>\n",
       "      <td>-0.006</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_7</th>\n",
       "      <td>0.74</td>\n",
       "      <td>0.96</td>\n",
       "      <td>-0.08</td>\n",
       "      <td>-0.59</td>\n",
       "      <td>0.044</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.1</td>\n",
       "      <td>0.042</td>\n",
       "      <td>0.0063</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_8</th>\n",
       "      <td>0.71</td>\n",
       "      <td>0.96</td>\n",
       "      <td>0.072</td>\n",
       "      <td>-0.53</td>\n",
       "      <td>-0.49</td>\n",
       "      <td>-0.22</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.42</td>\n",
       "      <td>0.14</td>\n",
       "      <td>0.015</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_9</th>\n",
       "      <td>0.7</td>\n",
       "      <td>0.96</td>\n",
       "      <td>0.092</td>\n",
       "      <td>-0.42</td>\n",
       "      <td>-0.53</td>\n",
       "      <td>-0.48</td>\n",
       "      <td>0.31</td>\n",
       "      <td>0.56</td>\n",
       "      <td>0.26</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.0036</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_10</th>\n",
       "      <td>0.66</td>\n",
       "      <td>0.92</td>\n",
       "      <td>-0.0096</td>\n",
       "      <td>0.093</td>\n",
       "      <td>0.29</td>\n",
       "      <td>-1.5</td>\n",
       "      <td>-1.5</td>\n",
       "      <td>0.44</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.73</td>\n",
       "      <td>0.18</td>\n",
       "      <td>0.016</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_11</th>\n",
       "      <td>0.63</td>\n",
       "      <td>0.92</td>\n",
       "      <td>-0.21</td>\n",
       "      <td>0.08</td>\n",
       "      <td>1.7</td>\n",
       "      <td>-0.69</td>\n",
       "      <td>-3.8</td>\n",
       "      <td>-1.8</td>\n",
       "      <td>1.7</td>\n",
       "      <td>2.1</td>\n",
       "      <td>0.92</td>\n",
       "      <td>0.18</td>\n",
       "      <td>0.014</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pow_12</th>\n",
       "      <td>0.62</td>\n",
       "      <td>0.9</td>\n",
       "      <td>-0.12</td>\n",
       "      <td>0.56</td>\n",
       "      <td>1.2</td>\n",
       "      <td>-2.9</td>\n",
       "      <td>-3.8</td>\n",
       "      <td>1.4</td>\n",
       "      <td>3.7</td>\n",
       "      <td>1.2</td>\n",
       "      <td>-0.53</td>\n",
       "      <td>-0.46</td>\n",
       "      <td>-0.11</td>\n",
       "      <td>-0.0096</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        rss intercept coef_x_1 coef_x_2 coef_x_3 coef_x_4 coef_x_5 coef_x_6  \\\n",
       "pow_1   4.1      0.75     0.54      NaN      NaN      NaN      NaN      NaN   \n",
       "pow_2   2.4      0.78     0.23    -0.14      NaN      NaN      NaN      NaN   \n",
       "pow_3  0.89      0.99     0.11    -0.56    -0.13      NaN      NaN      NaN   \n",
       "pow_4  0.85      0.96    0.043    -0.49   -0.038    0.021      NaN      NaN   \n",
       "pow_5  0.75      0.96    -0.13    -0.58     0.17     0.18     0.03      NaN   \n",
       "pow_6  0.74      0.95    -0.11     -0.5     0.17     0.11  -0.0093   -0.006   \n",
       "pow_7  0.74      0.96    -0.08    -0.59    0.044     0.16      0.1    0.042   \n",
       "pow_8  0.71      0.96    0.072    -0.53    -0.49    -0.22      0.4     0.42   \n",
       "pow_9   0.7      0.96    0.092    -0.42    -0.53    -0.48     0.31     0.56   \n",
       "pow_10 0.66      0.92  -0.0096    0.093     0.29     -1.5     -1.5     0.44   \n",
       "pow_11 0.63      0.92    -0.21     0.08      1.7    -0.69     -3.8     -1.8   \n",
       "pow_12 0.62       0.9    -0.12     0.56      1.2     -2.9     -3.8      1.4   \n",
       "\n",
       "       coef_x_7 coef_x_8 coef_x_9 coef_x_10 coef_x_11 coef_x_12  \n",
       "pow_1       NaN      NaN      NaN       NaN       NaN       NaN  \n",
       "pow_2       NaN      NaN      NaN       NaN       NaN       NaN  \n",
       "pow_3       NaN      NaN      NaN       NaN       NaN       NaN  \n",
       "pow_4       NaN      NaN      NaN       NaN       NaN       NaN  \n",
       "pow_5       NaN      NaN      NaN       NaN       NaN       NaN  \n",
       "pow_6       NaN      NaN      NaN       NaN       NaN       NaN  \n",
       "pow_7    0.0063      NaN      NaN       NaN       NaN       NaN  \n",
       "pow_8      0.14    0.015      NaN       NaN       NaN       NaN  \n",
       "pow_9      0.26     0.05   0.0036       NaN       NaN       NaN  \n",
       "pow_10      1.3     0.73     0.18     0.016       NaN       NaN  \n",
       "pow_11      1.7      2.1     0.92      0.18     0.014       NaN  \n",
       "pow_12      3.7      1.2    -0.53     -0.46     -0.11   -0.0096  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 系数矩阵为：\n",
    "coef_matrix_linear"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lasso回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "\n",
    "def lasso_regression(data, predictors, alpha, models_to_plot):\n",
    "    lassoreg = Lasso(alpha=alpha, normalize=True, max_iter=1e6)\n",
    "    lassoreg.fit(data[predictors], data['y'])\n",
    "    y_pred = lassoreg.predict(data[predictors])\n",
    "\n",
    "    # 绘制指定的alpha的图形\n",
    "    if alpha in models_to_plot:\n",
    "        myplot(data['x'], data['y'], y_pred, \\\n",
    "               models_to_plot[alpha], 'alpha=%.3g' % alpha)\n",
    "\n",
    "    # 记录模型拟合效果rss、截距和系数\n",
    "    return summary(data['y'], y_pred, lassoreg.intercept_, lassoreg.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 拟合了所有的 x\n",
    "predictors = ['x']\n",
    "predictors.extend(['x_%d' % i for i in range(2, pow_max)])\n",
    "\n",
    "# 设置正则系数\n",
    "alpha_lasso = [\n",
    "    1e-15, 1e-10, 1e-8, 1e-4, 1e-3, 1e-2, 1e-1, 1, 5, 10, 20, 50\n",
    "]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind = ['alpha_%.2g' % alpha_lasso[i] for i in range(0, len(alpha_lasso))]\n",
    "coef_matrix_lasso = pd.DataFrame(index=ind, columns=col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/chanson/anaconda3/envs/py37/lib/python3.7/site-packages/sklearn/linear_model/_coordinate_descent.py:476: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 0.35889189753831063, tolerance: 0.002422944771743981\n",
      "  positive)\n",
      "/Users/chanson/anaconda3/envs/py37/lib/python3.7/site-packages/sklearn/linear_model/_coordinate_descent.py:476: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 0.3588163321205684, tolerance: 0.002422944771743981\n",
      "  positive)\n",
      "/Users/chanson/anaconda3/envs/py37/lib/python3.7/site-packages/sklearn/linear_model/_coordinate_descent.py:476: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 0.35138644798246516, tolerance: 0.002422944771743981\n",
      "  positive)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEYCAYAAACQgLsAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9lUlEQVR4nO3deXhU9bnA8e87M1nYlwBhSwjIGhAUAq4oVVHAhavWDZVSrbT12vaq7S0WLfa6FOvWWtsqLgUUrFo3VNyoC6ggBGTf9wDKEnYSMpmZ3/1jJjAMsyaTmTMz7+d58pDMnDnnHT3vec9vOeeIMQallFIqnmzJDkAppVT60eKilFIq7rS4KKWUijstLkoppeJOi4tSSqm40+KilFIq7rS41JGIjBGRL+O9rFLpTnMnvWlxSTMiki0i/xaRzSJiRGRIHdd3h4iUikiViEwOeK/It43Dfj/31WV7SiWTiFwoIqtFpEJEPhORTmGWLfItU+H7zEUB798pIt+LyEEReVFEcvzee0BElomIS0Tur8evlDRaXNLTl8BNwPdxWNcO4EHgxTDLNDfGNPb9PBCHbSqVcCLSCngTuA9oCZQCr4b5yCvAt0AeMB74t4i09q3rEmAccCHQCegC/MHvs+uB/wXej++3sA4tLlESkXEiskFEDonIShG5MsRyRkR+KSIbRWSPiDwqIraAZR4TkX0isklEhvu9/mMRWeXbxkYR+WmscRpjnMaYPxtjvgTcQeLL8W1/q4jsFJFnRKRBmPW9aYx5GyiPNRalIHVyB7gKWGGMed0YcxS4H+gnIj2DxNod6A9MMMZUGmPeAJYBV/sW+RHwgjFmhTFmH/AAMKbm88aYKcaYD4BDtYgzJWhxid4GYDDQDO8ZyMsi0i7EslcCJXh3vpHALX7vnQGsAVoBfwJeEBHxvbcLuAxoCvwYeFJE+gOISKGI7A/zMyrK7zER6A6cBnQFOgC/j/KzoWwRkW0i8k/f2Z9S/lIld3oDS2o2Zow54ou9d5A4ewMbjTH+xWGJ37InrMv3e76I5IX43mlHi0uUfGczO4wxHmPMq8A6YFCIxR8xxuw1xmwF/gzc4PfeFmPMc8YYNzAFaAfk+7bxvjFmg/H6AvgYb1JijNlqjGke5md6pO/gS8SxwJ2++A4BDwPX1+a/CbAHGIi32T8AaAJMq+W6VJpKodxpDBwIiOcA3v06UKRlA9+v+T3YutKSFpcoichoEVlcc7YD9MF7BhVMmd/vW4D2fn8fGwcxxlT4fm3s28ZwEZknInt92xgRZhu10RpoCCz0+x4f+l5HRD7wG5i/MdLKjDGHjTGlxhiXMWYncAdwsYhkTAKpyFIodw7jbfn4a0rwrqtIywa+X/N72naDBdLiEgXfjJHn8B4884wxzYHlgIT4SIHf74V4B8UjbSMHeAN4DMj3bWNmzTZ8TfvDYX4iFgO8LY1KoLffWVszY0xjAGPMcL+B+dq0QGpusa37lQJSLndWAP381tsIOMX3eqAVQJeAE6l+fsuesC7f7zuNMRkzdqkHgeg0wnvg3A3ewUO8Z1+h/EZEWohIAfArws84qZEN5Pi24fINVl5c86avad84zM+xYuAbtM+tWa+I5IqIGGM8eBP9SRFp41u2g3hntgQlIg7fuuyA3bcuh++9M0Skh4jYfH3JTwGfG2MCuwtU5kql3HkL6CMiV/v2+d8DS40xqwM3aIxZCywGJvhy4kqgL94iBzAVuFVEikWkOXAvMLnm8yKS5duGDXD41mGP4rumDC0uUTDGrAQeB+YCO4FTga/CfOQdYCHene994IUotnEI+CXwGrAPGAXMqGXIa/C2UDoAH/l+r5mv/1u80yDnichBYBbQI8y67vV9fhze6c2VvtfAO73yQ7xN/eVAFSf2kasMl0q5Y4zZjXe210O+9ZyB33ikeGdWPuP3kevxTj7Yh3eizA9968AY8yHeSQefAVvxdvFN8Pvsc3hz6Qa805grgZtjjdnKRB8WFl8iYoBuxpj1yY5FqVSiuZNetOWilFIq7rS4KKWUijvtFlNKKRV32nJRSikVd45kbLRVq1amqKgoGZtWKqiFCxfuMca0TnYckWjuKCsJlzdJKS5FRUWUlpYmY9NKBSUiW5IdQzQ0d5SVhMsb7RZTSikVd1pcVHBl82HO495/lVLR09wBktQtpiyubD5MuQLcTrBnw49mQEGom9gqpY7R3DlGWy7qZJvneJPDuL3/bp6T7IiUSg2aO8docVEnKxrsPesSu/ffosGhl9UuAKWOizZ3MiBvtFtMnaxgkLc5v3mONzlCNeu1C0CpE0WTOxmSN1pcVHAFgyLv8MG6AOqSJGXzIxc0pawuUu5kSN7EpbiIyIt4n1+9yxgT7lkNKp3UdAHUnIH5dwHEusNnyNmcP82bDBUubyC23LFw3sSr5TIZeBrvA3JUpgjVBRDLDl+TSAe2xfdsLjVMRvMm84TrOos2d1Igb+JSXIwxs0WkKB7rUikmWBdAtM1+/0Sy2cHmAA+RJxGkCc2bDBaq6yya3EmRvEnYmIuIjAXGAhQWFiZqsyoeomimezyG8iNOjla7yWoxkDb2LMQN2LORUDu8fyJ5gAGjoVmB5fqOk01zJ4XF2j0cqcsMUiZvElZcjDGTgEkAJSUlep//VBHQTDej32FDbjFLyg6wbPsBVuw4wPZ9lew6VIXLc/x/a38Zx5m2VXxjerF16kE6NP+KUzs0o2/HZpxW0JyubRp7i45/IvUbdWLXmgUHKZNBcydF1WY8JJrZZuHypma7FsgdnS2mwts8B+N2IsaN21XFpClTeOTIpQA0yLLTu31TzuySR36zXNo2zaVhth0RQehHhdPF2Yeq6Hqwis3lR3jr2+28NM97n7uOLRpwUa98rrpwKn2ql2HrXMsxG6WsqrazwiLNNovHmE0CaHFRIS3ffoC5ZR24ydjJMoZqHBxseyYT+53KgE4t6NK6MXabRL0+j8ewcc8RFmzey6yVO3ll/lYmuzx0aH4611W34LpmR8lvmhv/qZpKJUM0XVy1VZcxmwSJ11TkV4AhQCsR2QZMMMa8EI91q8TyeAyfrNrJPz7fwOKy/WQ7mnGoy5P8V4uNdDx9KL8tOqvW67bZhK5tGtO1TWNuGFRIhdPFp6t38cr8rTzxyVr+8p91XNa3HXf3GkChf1I2yPNezZxmXWSaN2kuyouR9x1xsn73Ybbvq2THgUq+23+U/ZXVHDpazcHKaqrdBoO3N9Rhs9Ek10GjbAfNGmSR3yyXdr6fU1o3pkPhudgskjtJecxxSUmJ0WdSWIvbY3hn8Xb+9tl6Nuw+QmHLhtx6bmf+67QONGuYdXzBcP25dejr3bTnCNPmbWH6/K1UON3cfko5Yzpso02b9vDhuHpv5ovIQmNMSdxXHGeaOymsbD6V6z5nRXZfPj1cxNJtB1iz8xAFh5dxpm0V8zy9WGS607xhFi0aZtMk10GTXAfZdpuvqxmcbg+Hq1wcqXKxr6KaPYer8D+EN8iyc2mLMoY2XEer1u04fdVExF2N1FPuhMsb7RZTzFm3m4dnrmbVdwfp1a4pf73hdIb3aYvDHnDruXD9ubHOzw8oQJ1bNeLey4q5/Qdd+edXm5j8lYN/bMzj6YLPGeEb80l2M1+pWFW7PSzcso+1pZ9y7crbyTIueuPgEdd4KtsOYHTHnfx8yx+xearBkY1z1NvkdgnRO1A2HzZ/6cudc46tf9ehKrbvq2T9rsOs23WINd835+5thYwue4N+DicO8eB2VbFm7vu0zzuN5g2zE/LdtbhkghAH9K3lFUx5/TVyyr6msNFp/O7iTpzjmI+tZQ7Y25+8nnD9ubHOzw9RgFo2yubui3vwk3O78PfP1zPl6/VcYLeTLWCzZyFp2kWmLCpUazzwdb+/Xe1L+HpDOe8t3cFHK3ZyoLKaOxyfkOVwYRcPdnEzfaiLrCGDvfvyZhfgAXc1udu/hmDFJUTuZNltdGjegA7NGzCoc8tji7s9hm1Lc2DGO7g9Tpw4mLr4EK2W3sGR9mfRo+RCLu7dlpaN6q/QaHFJd0F2yur2Jbzw5SY+n/Ue/7Q9SE6WC3G/hXwl4HHB7BAtj3ADlLHOz4/QCmnWMIt7RvRi21mdePbtPJzrZ+O0tWDczN9iNzE28y0yNVOlmFAnQ4GvD5sIH47DuJ24xMFY+T2fHelM4xwHQ4vzuaR3WwY3aIH9lRngdiL2bLJOOc+7jWgH/WMcqLfbhE6nDYG8d2HzHLJzW/Kwr3vZuestRr31O+57pwdDi/O5pqSA87q1Dj45pw65o8Ul3QXslLuWzWLMm1Ws/O4gT7TbQu5+N2K8Z01eJvTOG26Asjbz86OYPdOxRUP+58c3saTsMpa9OgEOORHxeKdHR9NFZqGpmSrFhDqg+71u3E7KF7xOC1cVdjyIx3B5y41cd+UPGdKjNblZdt/K2gbPj2jvQF7bmWe+WWX2OY+DpxrwkCNu/n5uBc+ZIt76djszl31P+2a53HxWEaMGFR4fY61j7mhxSXe+ndK4nbjFwR1fN2JXbhXP3DSAYc1aw5Tpx28jga/lEm7nDTcHP9R7/mc/0SRSEP0KmtP32lG4J/8Ll8dJtcfOkytbcW2Pw3Rt0zj0By00NVOlmFAH9KLBGHs2xuXEaew8VtaD+7PnI7ixZ2Vx1VXXQ0Hbk9cXKj/C5VQccifwu4g9m7Z9h3JfQTG/HdaTWat28vK8LTzy4Wqe+s86fjigI2PP60JBHXNHZ4tlgANrv+STmW8wfWchzXucyxNnV9F85zfHk6Vmh/X/PV4H4Hi3HMrmU71hNjMOdGHCokZUVrsZ2a89Px9yCt3ym9R6+zpbTAXlf3AHXBtn83FFN95YuI0elUvY3Wog511wKcObb8Wx9SvL5064/F654yD//GoT7yzegccY/qfnPm7fehc2d3XI7YfLGy0uaW7ZtgP89KVS9hxxMn5EL0Z33IlMHZm4bqI5j8OnD3nPfsQOF4yHwXfHZdV7Dlfxj883MP2brVRWu7m4OJ/rBhZwTtdWft0RRNVvrMVFhePe8g1myhXgcVJtHDyYN5FLR4zkrFPyEIn+QuKY1GPuhPP9gaM888UGps/fSl+zhvHF5Zx+3uUxn5Rpt1ga+/fCbfzurWW0bpzDGz87m1M7NoM5byS2m6ger1Ju1TiH+y4r5r9/0JXJX29myteb+XjlThpk2Tm/e2t6tG1Cg2w7DbLakN/iRoYF66pQKoIv1+1h3Rsvc7PbO63XZnPz4Gn7ka6t6nfD9XmFfxhtm+Vy/xW9uX3IKTw7uxM5/TtC+6Yxr0eLSzrxnaF7Cs/l4WVNeP7LTZzVJY+nR51OXuMc7zKJ3mGjHbCsg5aNsrlraHfu6FbOtm8/4ZOKbkzemsWHK74/tsypHZoxrI8WFxVCkNbttn0V3D9jJbNW7eSSpr0Y7cjGeKq9V8B3TsCBPgG5A4Rs2bdpmst9lxXXerVaXNKFr3/WuJ24cLDo6D386Kyh3HdZ8YkXQyZqh/UXzSOT66psPtkvX0kXt5Of2rP56bCJeCrKcXY8myNt+uNOQvevShEBYxuum9/m+c2t+cusdQD8dlhPfnzOMOzflyR+Snt9506wadWV5XH5jlpc0oXf3YttxnBP8R4Gjgzx5NxEHOwTJdgT+VxVMPNubMaQa88mV6cfq3ACpha/9Mo0Ju4fxsXF+Uy4ojcdmjfwLpcheYMxcRmP1eKSJspbDaKRseMwBnFkM3DIyGSHVP9CPZFPBIzH+6PTj1UkvqnFHpcTp8fOl66ePDe6hKHF+cmOrH4kKG+0uKSBsr0VjHq3mkLzeyb2309B/0sy42Aa6ol8DfJOvNmlhR79qqxnTVYvnsn9P9rtKyW32/k8ce21J96sNd0kKG+0uKS4DbsPc9Pz31DhdPObn9xMQUFz7xuZcMuTcE/kyy9O/++v6sQYw0vztvDg+6tomtuZiTeO5KLifG/uLEzjfSdBeaPFJYVt2H2Y656dBxj+NfZMerXzTRfMlFueRLodTTp+ZxUX+yuc/Pr1pcxatZPzu7fmsWv60bpJTmbkToLyRotLitpSfoRRz83jVLOax0oOkedqBvh2iky65YkWERWjpdv28/TU6fSsXMwV5w7nshEDsdXctDFTcicBeaPFxUqi7Mratq+CUc99Q8/qVbxgexDb/GpY+JfjZ1lJuvhKqaSJIneMMUyfv5V33n2bqY4HyXG4kcUzoF+H45/R3IkbLS5WEelBXL7E2dWsL6Oe+4ZDR6t5ZOBBbKXVJ59lJeNaFqWSJYrccRaczb2lDXmtdBuP5m8m56A7+APoNHfiRouLVYRqjvsljrFn8cfcB9hzuJDpt51JW1rAt38Nfpal3UUqU0TIHeN2YoyD9VX38MsLRnB1j1bIS6+Ebp1o7sSFFherCNUc90scj8vQ4UApz46+ktMKmgN6lqVUuNypubDYbgwP9z9Az4t7eN/TvKl3WlysIlRz/NgFXlU4jYMzfjCSwd1an/g5TQ6VyULkTqn0prfHTpYYbI5sep454sTPaN7UKy0uVhJkhzcdB/J8lz+zb8Wn9D57BJdecGmSglPKwgJy56W5m5nwvpurWj/E/X330bjHD7SYJJgWF6vyDUS+e6ALDy1twtjz7uTSEb2SHZVSlubxGF56/XV2Lp3FT4rO4VdjbqFRjh7mkkH/q1uRbyDS46piqHFw+ymP8+thI058X/uLlTpBtdvD36ZO56eb7yQ7y4Vt9wxkV6+gsy41b+qfFhcr8g1E2vCQJS7u7Lbr+EVemXAFsVIxqnC6uH3aIorXzyYny4UNT8hZl5o3iWGLvIhKtD2tBlFl7LiwYXdkk3XKecffDDbtUqkMdqCimhuf/4bZa3dz2nmXYXPkeB8LHGLWpeZNYmjLxWIqnC7GzIKmnvv4y1lHaN3nwhPPsPQKYqWO2X2oiptf+IaNu4/w9xv7c3GfdlDcNuisS82bxNLiYiEej+HXry9hxY6DPD/6Olr3CvI8Cb2CWCkAvjtQyY3PfcOOA5W8MKbk+BT9YNOMNW8STouLhTz16TpmLvue343oyYXBCksNnaOvMtzW8gpGPT+PAxXVvHTrGQwsahn5Q5o3CaXFxSI+XP49f561jqv7d+S2wV2SHY5SlrV5j/eO4EecbqbddgZ9OzZPdkgqiLgM6IvIMBFZIyLrRWRcPNaZSVZ/f5C7XlvMaQXNeejKPohIskNSCaK5E5uNuw9z/aR5VFa7ma6FxdLqXFxExA78DRgOFAM3iEhxXdebKfYdcXLb1FIa5zh49uYB5GbZkx2SShDNndhs8BWWareHV8aeSe/2zZIdkgojHi2XQcB6Y8xGY4wT+BcwMg7rTXsut4dfvPItOw9U8ezNA8hvmpvskFRiae5EadOeI9wwaR4eY3hl7Jn0bNs02SGpCOIx5tIBKPP7extwRuBCIjIWGAtQWFgYh82mmCBXB//xg9V8uX4Pf/phX04vbJHkAFUSaO5EUjaffSs+5aHSprhMN/419ky65zdJdlQqCgkb0DfGTAImAZSUlJhEbdcSglwd/Obu9rzw5SbGnF3EtSUFyY5QWVjG5k7ZfDxTLqeJy8nTOPj+v16jSAtLyohHt9h2wP/o2NH3mqoRcHXwd0s+YdybyzirSx7jL9WbUWYwzZ0wDq76DI/LiQMPOeKm6NCiZIekYhCP4rIA6CYinUUkG7gemBGH9aaPmquDxY6xZzH524PclfsezwxxkbWjFOY87m3dqEyjuRPCroNHGbeoKdU4MGJH7NnQIO/EXCmbr7ljYXXuFjPGuETkDuAjwA68aIxZUefI0onv6mDXxtk8V3qA/zn4DLk2N/Lqa4CAx6U308tAmjvBlR+u4sbnv2F7RWe2XP4KPY8u8RaWD8cd71oeNvHEvzV3LCcuYy7GmJnAzHisK12ZjgMZ/00ueXufJifLhRgPuD017554B1eVMTR3TnSgopqbX5jP1r0VTLllED275AEXeVso/jeeXPXOyTei1NyxFL0rcoJM+Xozr5aW0bH/UL+7tmYd6y7Tm+mpTHe4ysWYyfNZt+sQk0aXcGaXvONv+nUtY8+GXiM1dyxOb/+SAF+u28MD769iaHE+1185AAZ1Oj4tGfRmeirjHa12c9uUUpZuO8DfRvXn/O6tT1wg2I0n84s1dyxMi0s9W7/rED+ftpBubRrz5HWneR/6FXgDPU0MlcGcLg+3T1vEvE3lPHFtP4b1aRt8wWB5o7ljWdotVo/2HnFyy+RSchx2nv9RCY31Wd5KncDtMdz52mI+Xb2LB0b24crTOyY7JBUnWlzqSZXLzU9fKmXnwaM8N3oAHVs0THZISlmKMYbfvbmM95d+x+9G9OSmMzslOyQVR1pc6oHHY3h66isMLJvMCxd69NYuSgUwxvDg+6t4tbSMX17QlbHnnXL8Tb1+JS1oP008lc3HbJrD++uruH3Lk+RkubB9NQO66hx8pfw9OWvdsdsf3Tm0+/F77wVez6LXr6QsLS7x4rt/mHFVMcwIdvFg0+tXlDrJpNkbeOo/67hmQEd+f1kxsm3B8XvviYDxeH80d1KaFpd42TwHj6sKGx6MCGKzgzE6B18pP9O+2cLDM1dzad92TLy6r3f2pP+994wNbDZANHdSnBaXOJld3ZOBxkG2uLA5cpBhE6GyXOfgK+Xz1rfbuPft5VzQsw1PXnsadpvvias1F0j639pFcyflaXGJVpDnsdR4d8kOfvUJ3NThEe7rsxf7KedpUigFx/JmrrsXd39kOLNzHn+/sT/ZDr+5RMEukFQpT4tLNII8j6UmAWYu+47/eXUxJZ1aMu6WS8jK1v+kSgHH8sbjquI04+C6thO590eXBH+Ut14QmXZ0KnI0Ap7HwuY5ALw8bwt3TF/E6QXNefHHA2mohUWp4zbPwfjGIbPExYS+e2mkFxJnDP0/HY2APmHT6Vye/GQtT/1nHRf0bMPTo07XwqJUgKWOU+lmHGSJC7sjG0fX85MdkkogPSJGw69P2FlwNvfNb8Crpeu4tqQjD195Kg67NgCV8jdvYzljZrq5pOmDPNz/AI26D9FurwyjxSVaBYMoa9SH/56+iKXbyvjFBV25a2h3RCTZkSllKfM37eWWyQvo2KIh9942hkZNcpIdkkoCLS5RmrVyJ3e9thgDPHvzAC7p7btza5hZZEplmq/X7+HWKaW0a57L9NvOoHW4wqK5k9bSu7hEs/OGWsb3+v78M/i/bxvz5rfb6d2+Kf+4cQCFeQ2PLxNiFplSKSvag37Acl+s3c3YqaUMa7aVh/rtp/H+htAkTN5p7qS19C0u0ey8oZYpm4/x3colBwdlrvHc8YPh3HFB1xOnUQabRaYJolJZtAf9gOUWnD+Z2z40XNqyjCcqf498VQ3zngz9ec2dtJe+I9Ehpg9HWqba7WHZV+8du5VLNi6eOaeSX1/S4+T5+YGPXtVbVahUF03eBCzncTn54qM36dmuCQ+fvh9xV0f+vOZO2kvflkvgLSUa5Hlv4+3f1PdbxmPP4osyN2v++As2V+Rwf7YDwY3NkU1enwtPXLd/d4BeWazSSTR547ecx+WkythplteW13vPJadJmxM/7180ArvbNHfSWmoWl2j6hP133iC38TYdB7IhpxerT3+Gg6s/Y2m5nQlrHuU8cUGDbGzD/4RU7g0+FhPYbTD47sR8b6XqKlLuRMibms94Ogzk5W5P8f3SWXRo34Gf7PsHMjvMvcFCdbdpUUlbqVdcYhkIrNl55zyOcTsRXxP+gxmv8eDBQ3x34CiQRdc2V/FA94/I3eJGjAc81VC5N3jR0L5ilaqizR2/vAm2r1c4Xdz56mI+WtGIm878Bb9u+THymd9yleUn547mTcZJveIS5U5a7fawcMs+Fmzay761efzGYycLQzV2/r23M/27tODcbq04t2srClo2hLIGMOXF4M15f4HdBtpXrFJFrAf4IPv6jv2V/GRKKau/P8i9l/bi1nM7I9sqYXaEnNC8yTipV1zC7KQVThefrNzJrFW7+HzNLg4ddSECPfK7MLX7U5zrWE2rPhfyYvHgky9+jLYPWPuKVaqK9QAfsK/PdZ7CL57+iqPVbl4YM5Af9GgTdLmI3W2aNxlBjDEJ32hJSYkpLS2t/QoC+o3X7jzEtHlbeHPRdg5VuWjVOJsLerbhgp75nNUlj2YNs8J+XikRWWiMKUl2HJHEO3ei4XJ7eOo/6/jrZ+u5tHkZE/ruo3WfCzV3VNi8Sb2WCxzrE166bT+PvzifL9buJttuY8SpbblhUCEDi1p6n3AXjF68pdLQkrL9uDweBnRqGX7BGAfRy/ZWcNdri1mweR939dzPL7bdj8x3wsK/aO6osFKnuPidca3N7sXjH6/hoxU7ad4wi99c0oPrBxaQ1ziKexjpwKJKQxM/WM2CzXsZf2kvxpxddLzbt5at9CqXm0lfbOTpz9bjsAlPXtePKw+/Cls0d1R0rFVcwt2KZcoVGLeTahzcU3UPa7OKufOi7txybhFNcrNCrzOQDiyqNPTcBR4+en8W095bw8It5/PI1X1ptGtRzK10j8cwa9VOJn64mo27jzDi1Lbcd1kx7Zo1gDLNHRU96xSXMN1VB1Z9RhPfFfM2Y7i96Hv63/jftGiUffI6dEBeZZqy+TT+11Vc5XZyRQMH1y+HEdsP8FSHz+jrm4IfqaXh3DyX1XM/4Lmy9ry7t4CivIZM/vFAhtQM2oPmjopJnYqLiFwD3A/0AgYZY2o/0uh7ap3gwe2qYt6st9lS3JaPV37PwXXZTMtykC0ubFnZXDj8aqgpLDUFJcwFXyfRi7dUksU7d/AVkSzgiUGHuX2rgz8sa8H0HDtZALYs9uYNopUxGAOHjro4tOErypd9ytK9dn64+2mKcfGoOLjm4smcff75wZ9TpLmjolTXlsty4Crg2TpHUjQYj6/JXY2Dx9e1ZtGaZeQ3zeHaISM42GEA+XsXhL7qVwSMx/uj/cHK+uKaO/7dVZ0GXMJ7Iwcyd2MvHv+kJVllXzGvqheLph4kN+tDqt2GfmYN07Ifpi0uikWwiwcbBru4OS9rNdgvrXNYKrPVqbgYY1YB8XlgVsEg7GPehc1zsBcNZlr+AMrXzKHd3lLsXZpCwXnAeSd+xn9w3tjAZgNE+4OV5cU7d07orgLkyyc4u2gwZ//sR+yvuIEhOw+zduchNu05Qm6WjSG7viJngwsbHgw2xGYHYxDNHRUnCRtzEZGxwFiAwsLC4Av5NbkblM2n4zvXe4vHnMeCd3MFDs4Hu6eRUikuptwJMnbZvGAQgzq3ZFBnv2nKZSNhywve7jTNHVUPIhYXEZkFtA3y1nhjzDvRbsgYMwmYBN4LwSJ+IJopwzrAqCwsKbkT7VR7zR1VzyIWF2PMRYkI5CTRThnWAUZlUUnJnVim2mvuqHpknanIgfTMSqnYad4oi6jrVOQrgb8CrYH3RWSxMeaSuEQGemal0la95o7mjbKApNy4UkR2A1ti/FgrYE89hFNbVosHrBeT1eKB0DF1Msa0TnQwsapF7qTS/4NksVo8YL2YYs6bpBSX2hCRUivdtdZq8YD1YrJaPGDNmOqTFb+v1WKyWjxgvZhqE0+QS3CVUkqputHiopRSKu5SqbhMSnYAAawWD1gvJqvFA9aMqT5Z8ftaLSarxQPWiynmeFJmzEUppVTqSKWWi1JKqRShxUUppVTcpUxxEZEHRGSpiCwWkY9FpL0FYnpURFb74npLRJonOZ5rRGSFiHhEJKnTGEVkmIisEZH1IjIuybG8KCK7RGR5MuNIFqvljtXyxheTJXLHSnnji6fWuZMyxQV41BjT1xhzGvAe8PskxwPwCdDHGNMXWAvck+R4ap4RMjuZQYiIHfgbMBwoBm4QkeIkhjQZGJbE7Seb1XLHankDFsgdC+YN1CF3Uqa4GGMO+v3ZCEj6TARjzMfGGJfvz3lAxyTHs8oYsyaZMfgMAtYbYzYaY5zAv4CRyQrGGDMb2Jus7Seb1XLHankDlskdS+UN1C13rHvjyiBE5CFgNHAA+EGSwwl0C/BqsoOwiA5Amd/f24AzkhSLwtK5o3lzXFrljaWKS6TnXxhjxgPjReQe4A5gQrJj8i0zHnAB06wQj8o8Vssdq+VNtDGp+LFUcYnh+RfTgJkkoLhEiklExgCXAReaBFw0lLTn68RmO1Dg93dH32uqnlgtd6yWN9HEZAFplTcpM+YiIt38/hwJrE5WLDVEZBjwv8AVxpiKZMdjIQuAbiLSWUSygeuBGUmOKWNZLXc0b0JKq7xJmSv0ReQNoAfgwXvL8Z8ZY5Ja1UVkPZADlPtemmeM+VkS4/F/Rsh+IL7P14ktlhHAnwE78KIx5qFkxOGL5RVgCN7bhu8EJhhjXkhWPIlmtdyxWt6AdXLHSnnji6fWuZMyxUUppVTqSJluMaWUUqlDi4tSSqm40+KilFIq7rS4KKWUijstLkoppeJOi4tSSqm40+KilFIq7rS4KKWUijstLkoppeJOi4tSSqm40+KilFIq7rS4KKWUijstLnEmImNE5Mt4L6tUJtE8Sn1aXDKAiFwoIqtFpEJEPhORTmGWLfItU+H7zEV+7/URkY9EZI+I6O20VUaJMY8eEJFlIuISkfsTGKZlaHFJcyLSCngTuA9oCZQS/pnlrwDfAnnAeODfItLa91418Bpwa70FrJQF1SKP1uN9INr79R+dNWlxqSURGSciG0TkkIis9D1sKNhyRkR+KSIbfWf8j4qILWCZx0Rkn4hsEpHhfq//WERW+baxUUR+WotQrwJWGGNeN8YcBe4H+olIzyCxdgf6430gUKUx5g1gGXA1gDFmje9BQStqEYdSJ0nHPAIwxkwxxnwAHKrFttKCFpfa2wAMBpoBfwBeFpF2IZa9EijBe+AeCdzi994ZwBq8T3r7E/CCiIjvvV14nzPeFPgx8KSI9AcQkUIR2R/mZ5RvHb2BJTUbM8Yc8cXeO0icvYGNxhj/hFgSYlml4iEd80ihxaXWfGcwO4wxHmPMq8A6YFCIxR8xxuw1xmzF+wjTG/ze22KMec4Y4wamAO2AfN823jfGbDBeXwAf401EjDFbjTHNw/xM962/MXAgIJ4DQJMgccayrFJ1lqZ5pNDiUmsiMlpEFtec4QB98J41BVPm9/sWoL3f39/X/GKMqfD92ti3jeEiMk9E9vq2MSLMNkI5jPeMzV9TgjfXY1lWqTpL0zxSaHGpFd8skeeAO4A8Y0xzYDkgIT5S4Pd7IbAjim3kAG8AjwH5vm3MrNmGrzl/OMzPjb5VrQD6+a23EXAKwcdNVgBdRMT/bKxfiGWVqpM0ziOFFpfaagQYYDd4BwzxnnGF8hsRaSEiBcCvCD/LpEY2kOPbhss3QHlxzZu+5nzjMD/TfIu+BfQRkatFJBf4PbDUGLM6cIPGmLXAYmCCiOT6Blf74k1OxCvXFxu+ZXKi+C5KBZOWeeT7Llm+5WyAw5cr9ijiTRtaXGrBGLMSeByYC+wETgW+CvORd4CFeA/c7wMvRLGNQ8Av8U793QeMAmbUItbdeGd7PeRbzxnA9TXvi8gzIvKM30euxztoug+YCPzQtw6ATkAlx8/WKvEOoioVszTPo+fw5scNeKf0VwI3x7rdVCbG6LVw9Um8Fxt2M8asT3YsSqUqzaPUoy0XpZRScafFRSmlVNxpt5hSSqm405aLUkqpuHMkY6OtWrUyRUVFydi0UkEtXLhwjzGmdeQlk0tzR1lJuLxJSnEpKiqitLQ0GZtWKigR2ZLsGKKhuaOsJFzeaLeYCq5sPsx53PuvUip6mjtAklouyuLK5sOUK8DtBHs2/GgGFIS6l6BS6hjNnWO05aJOtnmONzmM2/vv5jnJjkip1KC5c4wWF3WyosHesy6xe/8tGhx6We0CUOq4aHMnA/JGu8XUyQoGeZvzm+d4kyNUs167AJQ6UTS5kyF5o8VFBVcwKPIOH6wLoC5JUjY/ckFTyuoi5U6888ai4tItJiIvisguEVkej/WpFBGuCyDWZn/N2dynD3n/TePughqaNxkqUtdZbXLHgl1s8Wq5TAaeBqbGaX0qFYTqAoil2V/TWjmwLSPO5gJMRvMm84TrOou1y8zCXWxxKS7GmNkiUhSPdakUE6wLINpmv39i2Oxgc4CHyJMI0oTmTQYL1XUWS+5Y/KQsYWMuIjIWGAtQWFiYqM2qeIh1LKSm2V9zNhWqUPgnkgcYMBqaFeiYi8pc0eROipyUJay4GGMmAZMASkpK9FbMqaI2ze5oZ5sFJlK/USd2rengPqAnZikt1v04mtyJdFJmkdzR2WIqvNrObIlmtlk8xmwygJ6Ypaja7seRcifSSZlFckeLiwov2i6u2qrLmI1SVlZf+3G41o2FcideU5FfAeYCPURkm4jcGo/1Kguo2ZEvGJ+4s6DAqZoN8iw51bKuNG/SXCx3uohVwSAYfPfJ+Wih3EnKkyhLSkqM3jY8RYXrz41nX2/NuhrkwYfj6r2ZLyILjTElcV9xnGnupJhkjH8kMHfC5Y12i2U6/50fan/bimj7eqNNtprusjmPn9DMN4un49owG2fB2RzNH0C12+AxhppTJLsIDruQZbORk2Ujx2FDROr0n0ipWvPv9g2178e7AIXInUR3kWlxyQThdmr/KY0IeFyhi0O4/txo+nqDFCDTcSDlR5zs2F/Jjv2V7DpUxe5DVew5XEX5YSetDzTnPmPHgcFtBEpfwo4HGw5uc/6ORaZ72K9uE2iQZadfQXOm33Zm3f9bKgUn51SkAhHq5Ks2FxzX9pKAmi6yBLWitLiku3A77wkFweP7gAldHMIN7kcY+N9f4eTwtx/T3lWFDQ9uVxVTXp7KI0fKqXJ5TljWJpDXOIe8RtkcbFjMXzo8Tl/XMlq5dzFgzwxseLDhZnzvvazpdio2gZrGidsDLo+HarfhaLWbSqeblvsWU1w1B8psOjFAxSbYAT0wp4ZNjNz9FOrkqzYXHNfmkoAEdS/70+KS7sLtvP4FIbDlEmzwMdwsFd97ZtMcduUNpHRfActXrGbVdwexb19Aj8ol7DWNmZDlIAsXLnGwM28gN5/aiY4tGtC+ufcnv2kuLRtlY7f5d2Wd5f2nbD5M+QjcTmz2bAacfzkDCiJc91E2H6b8yvsdp7yQ8dOaVQxCHdADc2rVO5ELRKiTr9pccFybSwJq20VWhy47LS7pLtjO67/D+BcLiLwjBUwdrnS6WbJtP6Wb97JwCyzZ1pe9Rw4Bi3DYhMtbbuMR9x9wZFVj7NkcOP8hcjhEbufB3BPrQT7aizP9WWhqpkoxofadwJzqNRK2zA1fIELtu5H2af/B+bpcElCbSwrqeM2MFpd0F7jzwsk7zOC7T1w+jKPVbhZu2ce8jeXM3VDOkm37qXZ7h9O7tmnMhT3b0LegOf06NqNH2ybkzP0zfOoCPOCppqUcOnF7tfk+sRSH+r5OR6WvcGMWgQUhvzjmE7OIrwfrfqssr92YSRJOzLS4ZAL/nbcWzeONuw/z6epdzFm3h282lXO02oPdJvTp0Ixbz+3CoM4t6F/YguYNs0/+cLIP7rVJKqUg8phF4ElZvPetwIN7ZXlKnZhpcUkn0fSPRrHDuD2GBZv38snKnXy6eheb9hwBoEvrRlw/sJDB3VoxqHNLmuRmRY4pkQf3UN+/PhJfZYZkTutN8RMzLS7pItr+0RA7TLXbw1fr9/Dh8u/5ZOVOyo84yXbYOPuUPG45p4ghPdpQ0LJh7WJLxME9nl0ISgVKxoE+USdm4U5K65C7WlzSRSz9o74dxuMxfLOhnBlLdvDh8u/YV1FNkxwHF/Rqw7DebTmve2sa5Vh8Fwn2XAtXFcy8G4xJ+s37VIqI1OpPVvdqfZ+Y1eNJmcWPHCpqMZxZrd91iDcXbeftb7ez48BRGmbbGVqcz+V92zO4eytyHPYEBl4HoZ5rIQLG4/3RGWIqklha/emyHyXgpEyLS7oINivM72rcI1Uu3lu6g1cXlLFo637sNmFwt1b8dnhPhhbn0zA7BXeFUM+1CBx81RliKpxMm66eoJOyFDyiqJBqzqz8dh6PLYtnO/+Zp9e24IjTTdc2jRk/ohcjT29Pmya5yY64bsI91yKaqaFKQfIHzhMtQSdlWlxSXZC+YtfG2dh8t1nxuAwVaz5neN/buWFQAf0LW6TPjRwj3TFAi4qKRoRWf9pJ0EmZFhcrifVWCwF9xQeu/TdTyvJZ8nUDnjYOssQF9izG3jSaJt361X/8yaBFRMVDkFZ/2k4GSdBJmRYXq6jNTu3XvHW7qnj+pZf4a/UVnNd9ICt7vMzp7uXYOg+mSbolh1J1FepELlPGXxJwUqbFxSrC7dQhEmFDw9MpwIHNGKpx0LD7+cy65Dy6tmniW+KSxH8PpZIhllZ/uBO5TBt/qUdaXKwi1E4dJBHmu7ry9Gfrmb22knNz72Nspx2ces6l/LzHucn9DkolQ6yt/nAncnq7oLjR4mIVoXZqv0QwbievvjadcbuH0qpxNr8d1pObzrw4utuwKJWuYm31R2qd6DheXGhxsZJgO3XRYDy2LIzb4PTY+biiG/ddVsyoQYU0yE6Rix2Vqk8xtPqP5Zi2TuqdFherKpvPrmWzeH5rB0orxnFBg3V06n8xf7/kcnKztKgodUwUrf6g3V9aVOqVFhcL2rliNi3+/UNaeqq5UxzMPPNZhg//xfGr6OvwdDil0lKIVr8OziePFhcL2V/h5G+frSd33jR+ZavGIR7s4ubqlpvAv7Ck+zx8peJBu7+SSouLBVS7PUybt4UnZ63j4NFq7u51Hratb4O7Ggk848qUefhK1VZgyz7ClH5VP7S4JNnna3bxwHsr2bD7COd0zePeS4vp1a4plHULngja1FcqtFAte23xJ5wWlyTZUn6EB95byaxVu+jcqhHPjy7hwl5tjt/3K9zztrWpr1RwoVr22uJPOC0uCXa02s3fP1vPM7M34rAJ44b35JZzOpPtsEW/Ep3pogL84d0VrNxxMNlhJF03Z3PuMw4cGFzGwQOLm7Nu5dyQr6vIits3ZcLlvWP+XFyKi4gMA/4C2IHnjTET47HedPP5ml38/p0VbN1bwRX92vO7Eb1o2yzFb3uvlIWsyy7mgbw/0tu5lBXZfVmXXRz2dVV/6lxcRMQO/A0YCmwDFojIDGPMyrquO13sPHiU/3t3Je8v+44urRsx/SdncHbXVskOS6WR2pxZpq+zYnxd1Yd4tFwGAeuNMRsBRORfwEgg44uLx2N4ZcFWJn6wGqfLw68v7s5t53VJnccIK6VULcXQ0R9SB6DM7+9tvtcy2obdh7l+0jzGv7WcUzs04/PrGnJH1gxyvluY7NCUhYjIMBFZIyLrRWRcsuNRKl4SNqAvImOBsQCFhYWJ2mzCuT2G5+ds5IlP1pKbZedPP+zLNW12IFOv02mQ6gTapazSWTxaLtuBAr+/O/peO4ExZpIxpsQYU9K6des4bNZ61u86xNX/+Jo/frCa87u35pM7z+PakgJky5cnT4NUyq9L2RjjBGq6lJVKefFouSwAuolIZ7xF5XpgVBzWmzI8HsOLX23iTx+toVG2naduOJ3L+7Y7fs1K4IWPDfKOP6Mb9JqVzBWsS/mMJMViTXpVfcqqc3ExxrhE5A7gI7xTkV80xqyoc2QpomxvBXe/voT5m/YytDifh688ldZNck5cyP/CxwZ58OE4b6Gx2QEBj0u7y1RImdKlfBK9qj6lxWXMxRgzE5gZj3WlCmMM/164jftnrMAmwmPX9OPq/h2Ot1YC1Vz4OOdxvy4yT83a9KrhzBR1lzIwCaCkpMQkJjQL0KvqU5peoV8L+444Gf/2MmYu+54zOrfkietOo0PzBtF92L+LLLDlovcJyzQZ36Uclt5HL6VpcYnR3A3l3PnqYsqPVDFueE9uG9wFuy1EayWYwHuDgfYpZ6hM71KOKFiu1IxV1tyMUnPHsrS4RMnl9vCX/6zj6c/W0zmvEc//6Bz6dGhWu5UF3htMEyNjZWKXckxqciVw/GXYxONjlzoeY0laXKKwY38lv3zlW0q37OOaAR25/4reNMrR/3RKJUzg+Muqd3Q8xuL0CBnBp6t3ctdrS3C5DX+5/jRGnhblzQe0ya5U/ASOv/QaCVvm6niMhWlxCaHa7eGxj9fw7BcbKW7XlL/d2J/OrRqF/1BNQfGfbqxNdqXqLthzjPKL9QTOwrS4BLHr4FHumP4t8zfv5aYzC7n30mJysyLcbNK/T1gEjMf7o012peIj2Fil5pVlaXEJMG9jOXdM/5YjVa7YusH8+4SNDWw2QLTJrpTKSFpcfIwxvPDlJv74wWo6tWzI9NvOoHt+k+hXENgnPGwiVJZrk12pWOl4ZVrQ4gJUOF389o1lvLtkB8N6t+XRa/rSJDfrxIUi7fD6bHulak/HK9NOxheXreUVjH2plDU7D/G/w3rw8/NPOfkWLtHe40j7gJWKnY5XpqWMLi5fr9/D7dMXYQz8c8xAhvRoE3xBvceRUvVHxyvTUkYWF2MMU+du4f/eW0mXVo14bnQJReGmGes9jpSqPzpemZYyrrg4XR4mzFjOK/PLuKhXG5687rSTx1cC6XiKUvVH8ystZVRx2V/h5GcvL2Texr3cPuQUfn1xD2zR3nQy1HiKzmxRKrRo80PHK9NOehcXvx17Y24xt04pZfu+Sp64th9X9e940jIx79z6MCOlQtP8yGjpW1z8dmy3LYsJ7ns5YO/J9NvOoKSo5UnLnLTzR1N0dKBfqdA0PzJa+hYXvx3buAxDctfw8E9vpaBlw6DLnLDzR3vGpQP9SoWm+ZHR0ra4mE7n4hIH4jG4JYsbzj+dhiv+cWJLJNTOH+mMy79VowORSgWnA/UZLTWLS4Quq2q3h/Hzc1lfeQ+3FuzgkpJe5Hz8u5NbIqGedNcgL/QZV7BWzeC7E/TFlUox4QbqdTJMWku94hKhy+pIlYvbpy3ii7W7+dWFlzLiom7Il0+EbomEe9JdsLn22o+s0l20B32dDKPCSL3iEubgvvtQFbdMXsDK7w4y8apTuX5Qofcz0fT9Bq63sjx4i0T7kVU6i/agr5NhVASpV1xCHNy3lB9h9Ivz2XWwiudGD+CCnvnHPxNN32+0RUP7kVU6i/agr5NhVASpV1yCHNyXbz/AmH8uwO3xMP22Mzi9sEXwz4U7s4qlaOgFXypdBR70G+R5xyEDc6I2k2EC805P0tJa6hUXOOHgPndDObdNLaVproOpY8+ma5vG4T8b7sxKi4bKdP4H/XC3v491MkyovNN8S1upU1yCtDY+WbmT/56+iE4tGzL11kG0a9Yg8nq0r1dlsmjGQ2oO+nMeD58rsUyG0bzLONYqLqF2/CBnPW/ubs9v/r2UPh2aMXnMQFo0yo5uG9rXqzJVrDO0os2VaCbDaN5lHOsUl3A7fsDOu/CLd7lr+ZmcfUoek0aX0DjHcXwd0ZyVaV+vykSxth6izZVoCofmXcapU3ERkWuA+4FewCBjTGmtVxZux/fbeavFwUMrWjK0OJ+/3nA6ud8vjP3xqNrXqzJRbVoPgblSl8kwmncZpa4tl+XAVcCzdY4k3I5fMAgz+h3mzHqbP69rQ0G/ITx2TT+ydpTq41GVilZdWw86GUbFoE7FxRizCjj5mfO1EWzH950lmU7n8sCSJry49lyuH1jAQ1eeit0m+nhUpSIJbGnUpQjooLyKQcLGXERkLDAWoLCwMPhCgdeiTLkC43ZSjYPFR+9hzNkXM+Hy4uPFTB+PqlRo8b7Fig7KqxhELC4iMgtoG+St8caYd6LdkDFmEjAJoKSkxET8wOY5GLcTMW5sxvDLU3Zyvn9hAR0kVKpGsLGQeLc0NN9UDCIWF2PMRYkIJJC78FzcOLAZg7Fncf7FVwbvftO+XpXpQrVQYmlp6OOIVZxZZyqyH7fH8Jt5OWw+eg93d9/NORf9l+7QKq0kZKZlpJZGTUGJZaalUlGq61TkK4G/Aq2B90VksTHmkrqs0+0x/Ob1Jbz57XbuHnoZ51zYrS6rU8qq6m+mZeD9wCLd1VhnWqp6UNfZYm8Bb8UploDC0p1faGFRaareZlpG2wrRmZaqnlmqW+yNhdu0sCgVIKaZlpHuB1ZDZ1qqemap4vLDAR3Ja5zNhb3yIy+slMUlZaalPpdIWYSliovNJlpYVNpIykxLfS6RsghLFRelVBxo0VAWYEt2AEplIhG5UkS2AWfhnWn5UbJjUiqexJjIF8vHfaMiu4EtMX6sFbCnHsKpLavFA9aLyWrxQOiYOhljWic6mFjVIndS6f9BslgtHrBeTDHnTVKKS22ISKkxpiTZcdSwWjxgvZisFg9YM6b6ZMXva7WYrBYPWC+m2sSj3WJKKaXiTouLUkqpuEul4jIp2QEEsFo8YL2YrBYPWDOm+mTF72u1mKwWD1gvppjjSZkxF6WUUqkjlVouSimlUoQWF6WUUnGXMsVFRB4QkaUislhEPhaR9haI6VERWe2L6y0RaZ7keK4RkRUi4hGRpE5jFJFhIrJGRNaLyLgkx/KiiOwSkeXJjCNZrJY7VssbX0yWyB0r5Y0vnlrnTsoUF+BRY0xfY8xpwHvA75McD8AnQB9jTF9gLXBPkuOpeUbI7GQGISJ24G/AcKAYuEFEipMY0mRgWBK3n2xWyx2r5Q1YIHcsmDdQh9xJmeJijDno92cjIOkzEYwxHxtjXL4/5wEdkxzPKmPMmmTG4DMIWG+M2WiMcQL/AkYmKxhjzGxgb7K2n2xWyx2r5Q1YJncslTdQt9xJqRtXishDwGjgAPCDJIcT6Bbg1WQHYREdgDK/v7cBZyQpFoWlc0fz5ri0yhtLFZdIz78wxowHxovIPcAdwIRkx+RbZjzgAqZZIR6VeayWO1bLm2hjUvFjqeISw/MvpgEzSUBxiRSTiIwBLgMuNAm4aCgpzwiJ3XagwO/vjr7XVD2xWu5YLW+iickC0ipvUmbMRUT8n3s8ElidrFhqiMgw4H+BK4wxFcmOx0IWAN1EpLOIZAPXAzOSHFPGslruaN6ElFZ5kzJX6IvIG0APwIP3luM/M8YktaqLyHogByj3vTTPGPOzJMZzJfBXoDWwH1hsjLkkSbGMAP4M2IEXjTEPJSMOXyyvAEPw3jZ8JzDBGPNCsuJJNKvljtXyBqyTO1bKG188tc6dlCkuSimlUkfKdIsppZRKHVpclFJKxZ0WF6WUUnGnxUUppVTcaXFRSikVd1pclFJKxZ0WF6WUUnH3/4Afexed3/+PAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "models_to_plot = {1e-15: 221, 1e-3: 222, 1e-2: 223, 1e-1: 224}\n",
    "\n",
    "for i in range(len(alpha_lasso)):\n",
    "    coef_matrix_lasso.iloc[i, ] = lasso_regression(data, predictors,\n",
    "                                                   alpha_lasso[i],\n",
    "                                                   models_to_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rss</th>\n",
       "      <th>intercept</th>\n",
       "      <th>coef_x_1</th>\n",
       "      <th>coef_x_2</th>\n",
       "      <th>coef_x_3</th>\n",
       "      <th>coef_x_4</th>\n",
       "      <th>coef_x_5</th>\n",
       "      <th>coef_x_6</th>\n",
       "      <th>coef_x_7</th>\n",
       "      <th>coef_x_8</th>\n",
       "      <th>coef_x_9</th>\n",
       "      <th>coef_x_10</th>\n",
       "      <th>coef_x_11</th>\n",
       "      <th>coef_x_12</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha_1e-15</th>\n",
       "      <td>0.72</td>\n",
       "      <td>0.96</td>\n",
       "      <td>-0.043</td>\n",
       "      <td>-0.55</td>\n",
       "      <td>-0.052</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.062</td>\n",
       "      <td>0.0016</td>\n",
       "      <td>0.00032</td>\n",
       "      <td>0.00061</td>\n",
       "      <td>-0.00014</td>\n",
       "      <td>1.1e-05</td>\n",
       "      <td>1.7e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_1e-10</th>\n",
       "      <td>0.72</td>\n",
       "      <td>0.96</td>\n",
       "      <td>-0.043</td>\n",
       "      <td>-0.55</td>\n",
       "      <td>-0.052</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.062</td>\n",
       "      <td>0.0016</td>\n",
       "      <td>0.00032</td>\n",
       "      <td>0.00061</td>\n",
       "      <td>-0.00014</td>\n",
       "      <td>1.1e-05</td>\n",
       "      <td>1.7e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_1e-08</th>\n",
       "      <td>0.72</td>\n",
       "      <td>0.96</td>\n",
       "      <td>-0.044</td>\n",
       "      <td>-0.55</td>\n",
       "      <td>-0.05</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.062</td>\n",
       "      <td>0.0015</td>\n",
       "      <td>0.00027</td>\n",
       "      <td>0.0006</td>\n",
       "      <td>-0.00014</td>\n",
       "      <td>1e-05</td>\n",
       "      <td>1.7e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.0001</th>\n",
       "      <td>0.78</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.017</td>\n",
       "      <td>-0.46</td>\n",
       "      <td>0</td>\n",
       "      <td>0.012</td>\n",
       "      <td>-0.0078</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-9.3e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.001</th>\n",
       "      <td>1.1</td>\n",
       "      <td>0.87</td>\n",
       "      <td>0.082</td>\n",
       "      <td>-0.3</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0.0015</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>-0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.01</th>\n",
       "      <td>2.7</td>\n",
       "      <td>0.72</td>\n",
       "      <td>0.2</td>\n",
       "      <td>-0.13</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.1</th>\n",
       "      <td>24</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_1</th>\n",
       "      <td>24</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_5</th>\n",
       "      <td>24</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_10</th>\n",
       "      <td>24</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_20</th>\n",
       "      <td>24</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_50</th>\n",
       "      <td>24</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "      <td>0</td>\n",
       "      <td>-0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rss intercept coef_x_1 coef_x_2 coef_x_3 coef_x_4 coef_x_5  \\\n",
       "alpha_1e-15  0.72      0.96   -0.043    -0.55   -0.052     0.11     0.16   \n",
       "alpha_1e-10  0.72      0.96   -0.043    -0.55   -0.052     0.11     0.16   \n",
       "alpha_1e-08  0.72      0.96   -0.044    -0.55    -0.05     0.11     0.16   \n",
       "alpha_0.0001 0.78      0.95    0.017    -0.46        0    0.012  -0.0078   \n",
       "alpha_0.001   1.1      0.87    0.082     -0.3       -0        0       -0   \n",
       "alpha_0.01    2.7      0.72      0.2    -0.13        0       -0        0   \n",
       "alpha_0.1      24      0.16        0       -0        0       -0        0   \n",
       "alpha_1        24      0.16        0       -0        0       -0        0   \n",
       "alpha_5        24      0.16        0       -0        0       -0        0   \n",
       "alpha_10       24      0.16        0       -0        0       -0        0   \n",
       "alpha_20       24      0.16        0       -0        0       -0        0   \n",
       "alpha_50       24      0.16        0       -0        0       -0        0   \n",
       "\n",
       "             coef_x_6 coef_x_7 coef_x_8 coef_x_9 coef_x_10 coef_x_11 coef_x_12  \n",
       "alpha_1e-15     0.062   0.0016  0.00032  0.00061  -0.00014   1.1e-05   1.7e-05  \n",
       "alpha_1e-10     0.062   0.0016  0.00032  0.00061  -0.00014   1.1e-05   1.7e-05  \n",
       "alpha_1e-08     0.062   0.0015  0.00027   0.0006  -0.00014     1e-05   1.7e-05  \n",
       "alpha_0.0001        0       -0        0        0        -0         0  -9.3e-07  \n",
       "alpha_0.001    0.0015       -0        0       -0         0        -0        -0  \n",
       "alpha_0.01          0       -0        0       -0         0        -0         0  \n",
       "alpha_0.1          -0        0       -0        0        -0         0        -0  \n",
       "alpha_1            -0        0       -0        0        -0         0        -0  \n",
       "alpha_5            -0        0       -0        0        -0         0        -0  \n",
       "alpha_10           -0        0       -0        0        -0         0        -0  \n",
       "alpha_20           -0        0       -0        0        -0         0        -0  \n",
       "alpha_50           -0        0       -0        0        -0         0        -0  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.options.display.float_format = '{:,.2g}'.format\n",
    "\n",
    "coef_matrix_lasso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ridge回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "\n",
    "def ridge_regression(data, predictors, alpha, models_to_plot):\n",
    "    ridgereg = Ridge(alpha=alpha, normalize=True)\n",
    "    ridgereg.fit(data[predictors], data['y'])\n",
    "    y_pred = ridgereg.predict(data[predictors])\n",
    "\n",
    "    # 绘制指定的alpha的图形\n",
    "    if alpha in models_to_plot:\n",
    "        myplot(data['x'], data['y'], y_pred, \\\n",
    "               models_to_plot[alpha], 'alpha=%.3g' % alpha)\n",
    "\n",
    "    # 记录模型拟合效果rss、截距和系数\n",
    "    return summary(data['y'], y_pred, ridgereg.intercept_, ridgereg.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 拟合了所有的 x\n",
    "predictors = ['x']\n",
    "predictors.extend(['x_%d' % i for i in range(2, pow_max)])\n",
    "\n",
    "# 设置正则系数\n",
    "alpha_ridge = [1e-15, 1e-10, 1e-8, 1e-4, 1e-3, 1e-2, 1e-1, 1, 5, 10, 20, 50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEYCAYAAACQgLsAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABBYElEQVR4nO3dd3xb5dXA8d+RZHkmcWJnEjvOcCYkITEOK5QwwwwpUGbZpNDy0gHlpYUCfSktLaSlg9KGBhL2hgAJoWXGoUAGJCGD7OHs7REPWdLz/iE5URRJlm3ZupLO9/PRxxp3PEru0bnPuM8VYwxKKaVULNniXQCllFLJR5OLUkqpmNPkopRSKuY0uSillIo5TS5KKaViTpOLUkqpmNPk0koicp2IzI31skolO42d5KbJJcmIiFNEXhORDSJiROTUVm7vNhFZICL1IjIt6LMi/z6qAx6/as3+lIonETldRL4VkRoR+VhE+kRYtsi/TI1/nTOCPv+piGwXkUoReUpE0gM+e1BEvhERt4g80IZfKW40uSSnucDVwPYYbGsr8BvgqQjL5BpjcvyPB2OwT6XanYjkA28AvwK6AAuAlyOs8iLwNZAH3AO8JiJd/ds6G7gbOB3oA/QDfh2w7hrgLmBmbL+FdWhyiZKI3C0ia0WkSkSWi8jEMMsZEbldRNaJyG4ReUREbEHLPCoi+0RkvYicE/D+9SKywr+PdSLyg+aW0xjjMsY8ZoyZC3hClC/dv/9NIrJDRP4hIpkRtveGMeYtYE9zy6IUJE7sAN8FlhljXjXG1AEPACNEZHCIsg4ERgH3G2NqjTGvA98AF/sXuRaYaoxZZozZBzwIXNe4vjFmujHmPaCqBeVMCJpcorcWGAt0wncG8pyI9Ayz7ESgBN/BNwG4IeCzMcBKIB/4AzBVRMT/2U7gfKAjcD3wJxEZBSAihSKyP8Ljyii/x8PAQGAkMAA4CrgvynXD2Sgim0Xkaf/Zn1KBEiV2hgGLG3dmjDngL/uwEOUcBqwzxgQmh8UByx62Lf/z7iKSF+Z7Jx1NLlHyn81sNcZ4jTEvA6uB0jCL/94Ys9cYswl4DLgi4LONxpgnjTEeYDrQE+ju38dMY8xa4/Mp8G98QYkxZpMxJjfC44WmvoM/ECcBP/WXrwr4LXB5S/5NgN3Acfiq/aOBDsDzLdyWSlIJFDs5QEVQeSrwHdfBmlo2+PPG56G2lZQ0uURJRK4RkUWNZzvA0fjOoEIpD3i+EegV8PpgP4gxpsb/NMe/j3NE5AsR2evfx7kR9tESXYEsYGHA95jtfx8ReS+gY/6qpjZmjKk2xiwwxriNMTuA24CzRCRlAkg1LYFipxpfzSdQR0I3XTW1bPDnjc+TthksmCaXKPhHjDyJ78czzxiTCywFJMwqBQHPC/F1ije1j3TgdeBRoLt/H7Ma9+Gv2ldHeDSZDPDVNGqBYQFnbZ2MMTkAxphzAjrmW1IDaZxiW48rBSRc7CwDRgRsNxvo738/2DKgX9CJ1IiAZQ/blv/5DmNMyvRd6o9AdLLx/XDuAl/nIb6zr3B+LiKdRaQA+DGRR5w0cgLp/n24/Z2VZzV+6K/a50R4HEwG/k77jMbtikiGiIgxxosv0P8kIt38yx4lvpEtIYmIw78tO2D3b8vh/2yMiAwSEZu/LfkvwCfGmODmApW6Eil23gSOFpGL/cf8fcASY8y3wTs0xqwCFgH3+2NiIjAcX5IDeAa4UUSGikgucC8wrXF9EUnz78MGOPzbsEfxXROGJpcoGGOWA5OBz4EdwDHAZxFWmQEsxHfwzQSmRrGPKuB24BVgH3Al8HYLi7wSXw3lKOB9//PG8fr/i28Y5BciUgl8AAyKsK17/evfjW94c63/PfANr5yNr6q/FKjn8DZyleISKXaMMbvwjfZ6yL+dMQT0R4pvZOU/Ala5HN/gg334Bspc4t8GxpjZ+AYdfAxswtfEd3/Auk/ii6Ur8A1jrgW+39wyW5nozcJiS0QMUGyMWRPvsiiVSDR2kovWXJRSSsWcJhellFIxp81iSimlYk5rLkoppWLOEY+d5ufnm6KionjsWqmQFi5cuNsY0zXe5WiKxo6ykkhxE5fkUlRUxIIFC+Kxa6VCEpGN8S5DNDR2lJVEihttFlNKKRVzmlxUaOXzoGyy769SKnoaO0CcmsWUxZXPg+kXgscFdidc+zYUhJvEVil1kMbOQVpzUUfaUOYLDuPx/d1QFu8SKZUYNHYO0uSijlQ01nfWJXbf36Kx4ZfVJgClDok2dlIgbrRZTB2poNRXnd9Q5guOcNV6bQJQ6nDRxE6KxI0mFxVaQWnTB3yoJoDWBEn5vKYTmlJW11TspEjcxCS5iMhT+O5fvdMYE+leDSqZNDYBNJ6BBTYBNPeAT5GzuUAaNykqUtxA82LHwnETq5rLNOBv+G6Qo1JFuCaA5hzwjYFUsTm2Z3OJYRoaN6knUtNZtLGTAHETk+RijJkjIkWx2JZKMKGaAKKt9gcGks0ONgd4aXoQQZLQuElh4ZrOoomdBImbdutzEZFJwCSAwsLC9tqtioXmNnE1Ve1vFBhIXmD0NdCpwHJtx/GmsZPA2iJ2EiRu2i25GGOmAFMASkpKdJ7/BLH32zI6vXoJ4nHhkTSmFf+FDZnD6NExg4IuWRR0yaS4ewc6ZqQdWina0WbBgTTiysOb1izYSRkPGjsJKsomrroGD9sr6qiqc1Pl6kfaKdPosnseB3oeT71nAJ13VtGjUyY56f6f60hx07hfC8SOjhZTR9h7wMVTc9cz85ttnLPvBX7mcGEXL15vA7WrP2WWdGFfTcNh6/TNz2ZYr44M6dmRgd07MLD7MHqfdBx2mxyx/Xq3h/K9tWys7kPtsf+gw/YvWOI4hq8+gI6ZXzMmbQ3fW3YbdtOAWKyTUqmoBTVxmfVlrMsYyuLy/Swu38+KbVVs2lvD9sq6ECuPBhqAzw++k5PuoKBLFkN6duDUkU8w0rOU7sPPIL0lfTbtQJOLOmhnVR1PzlnHc19sos7t4TsDuzJo0DnYFs3AeBtwOJzcfu313F5QSq3Lw+Z9NWzcU8O32yv5ZksFX2/az7tLth22zY4ZDnKznDgdNqrr3FTVNXDA5QlYIo00+yl0zUmnS049q3dW07vqP2B3IeLF467HteoTMjW5qERTNBZjd2LcLtw4uPXTDD6c9SkAWU47Q3t25OTifAq7ZHFUbiadMtPITneQ5bTjNQaX20u928u+GhfbK+rYVlHH+t0HmLt6N29UOYFRpM/fz3FFX3LigDzGD+tBv1gPc26FWA1FfhE4FcgXkc3A/caYqbHYtmof7yzeyi/f/IYal4cJI3rxw3H9GdCtg+/DkUcdUc3OdNop7t6B4u4dOGNo94PbqaprYPXOalZtr2JrRR2VtQ3sr3FR7/bSIcNBh4w0OmWmUdgli8K8LAq7ZJGX7UTkUA3HvTEL27Mz8LhduIyD35ftYsL2X3HsKRckVQ1G4yY5eb2GL9bv4bWFTrbV/5JjvUtZ4jiGTgNO4OHifI4t7MyAbjkha/XR2lNdz9eb9vPftXv479rd/GH2Sv4weyUTu3bgD+LAAb5af2aebyaAODSRxeU2xyUlJUbvSWENVXUN3P/2Mt74agvHFuYy+dIR9OuaE36FSO25sWzr9W9re0M2ncvuw+5tAHsajuvfbZMgEZGFxpiSmG84xjR2rOtAvZtXFpTz1GfrKd9bS4d0B+eP6MUFI3pyXFEX0rYuCD/8uJVxs62illnfbGfWN9swm77kBPsKuvfoxVX7nsDmbbvm5Uhxo81iKWztrmpumDaf8r01/Pj0Yv7ntAE47BGmm4vUntvc8flNBZJ/qGaPsskY3Ih4cXsaWPLZuwy/PHlqLyrx7a9x8WSZrzm5oraBkj6dufOsQZw9rAcZaXbfQuHioyXXhIWInZ6dMrnx5L7ceHJfyveO5IV5m6j88s8Yr6952et2YdbNwd6OtRdNLqkgxEG5cOM+bpw+n5Gs4vkT9tN7cFcId2bVKFJ7bnPH50d7JlU0FrE7MR4XHnHw4jfVZDgeYOCYc5OqiUxZVLgf9PJ5NKydw4z9/fj1omyKXct5oEc5Q849l8HHnXjkdsLFR0uuCWsidgq6ZPG/4wfTMPhaeOYN3F4XDcbOX+bu4pyt9zHsxPOw9xkTo3+g8DS5JLsQB+X7lYXc/uLXnJ6zgb+5H8T2dQMs/jMg4HWHP3gjjcFv7vj8aDsb/cOaZUMZOHN54L27sX/jxrviCWzXvRNVgtn2zad03jWPjOLvaEJS0Qvzg242fYln2oWI18V5xsGBvB/x/YonsO1rgPdfgB7NiJ2WXBMWZeykFR0P17+Dd30Z6yvS+MnCB7F/66Zh5RMsOuM5Rp101mF9nWH/DVrYZKfJJdkFHZRLP5vJrYtLGd47l0cGV2Gb2+D/zOtfwYQ/eCNdvxLNtS3RBlIwfxNZetlkjHgQvLg9Lrzr5uCIYu6lLm9cisM0YP47GdFhzSpaIX7Q16QPZcGrL3KJx4VDvNhtHq7tvBj2NUT+4Q8XHy29JqwZsWMrKGVoQOyAmw/fe53HVuby6wuHhe9jbeWwZk0uyS7goHRLGvctyeWE/nk8eU0JWTvS4PM/HppGIrDmEu7gjTTja7jPAs9+ogmkCN9F7E68bl81/83dfbiyiVWqVn5MprcBu3jjPjRTJZiA2DF2J6/sKuLe2XMYk9afS+xOTON1WEMmwMbPm/7hDxcfkWIqxrGDx4Xd7mTYcefy3Pz9jH+sjFtP7c+tp/Y/1D/UqJXDmjW5JDv/mdE3n73L/Us6k9XvRKae5iX9yz8fecBC7K/sDXX2M/aOVn0X24Yypqztzl8WZjO0dD8jC3LDrvJhTTFn48AuHl9wWWjuJWVx/uNtz7IP+f2KfNbM38Qfe29h7JkX4ch85/BY6T40IWKHDWVI0VjOKyil9JR6Hpq5nD9/uJq3F2/l0UuHM7pPl0PrtLS25KdDkVPA6ws3c+drizmpfz5TT/eS/vzE9ruCt2wyfPSQ7+xH7HDaPS0PkAAVtQ2Mf2wOWU47M28fe+RZF2CM4ezH5jDKtpqHR1VEDHwdiqyCGWN47stNPDRzOaWONTxt+41/SHw7XfneRrETbO7q3dz9xhK27q/llu/058dnFJPuCBjlFiFpRoobvc1xkntn8VZ+/tpiTuyfx7+uLSF983/b9x7fzbllcjN0ykzjD5cMZ+2uAzz47vKQyyzbWsmqHdUcPeZMX1Bqc5iKUnW9m9tfWsSv3lrKmL55/P2kGl9iaa+4gTaLnWAnF+cz+yen8L2SAv7+yVom/O0zVu+o8n1YUNri2NFmsWQSdJbx/rLt/OTlRZT06cKT15T4zu5bWdVttmg7LFtgbHFXJp3Sjylz1nHSgHzOzS0/bD+vLdyM027jguG9YrZPlaQCYmeVcwi3PreQ9bsPcNf4QdxySn9sW4Av/tR+cQNtGjuHKZ9HzoYyHi4dy5lDS/jf15dw4d8+4zcXHc3Fo3u3eLOaXJJFUPvswnHTuW2Wh+G9O/HU9ceR5fT/V7fXARsoUodlK9151iC+XL+XF15/jfGO32Dz+Jot3Gf9lryv53Nz3xPolJXW9IZU6gqIHY8tjfsb7qHCOZTnbzqeE/rn+ZaJR9w07rct9xX0u3H6+If5+Pgd/H5FPne86uHL9Xv49YVHk+k8stm5KZpckkXAyA6vx8Uns99gUI+rmXZ96aGpuhu19QHbjpzbFjB9QBnv71yAcbsAL7jrsb33c271epFtb0D54KT5vqoNbCjDeFyI8WDchnNy1vDYLZPo3jHj8OWSKG5C3snSXQ+z7qCDMTxodzL0uL9yz4LN9O6cxe2nFzd7F5pckoW/ucvrcVHvtbMldzTP3TiGTplJfNbuP+vK9bi4xG7H7bFhACOCzXhwiMF4G3T4sYqooeAkDA5sxuC1pXHZpVeQHpxYkkm4O1mKgPGC8SIeF1d220T/m69jRITRmJFockkWBaUsOeMZPpj1OuUdR/OrW64jN8sZ71K1rYDamh3YUfw91rhyWVudzhV7/47gwabDj1UE1fVubv1QOFD3C+4YuJMTT78IKWz7qVHiKtydLDPzYPbdh/UrjSnIa/FuNLkkiU9X7WLSux765F3JizcfT5e9i+Cr+N+Nrk0FDU7o9Z0b6FVQyikA5edZ4m58yrp2VdVzw7T5LN9WycMXX8JJJQW+DyxyJ8c2E+lOljG8XkeTSxL497Lt3PbC1/TvlsNzN5aSt2+xZe5G16aamo4mGb+zionyvTVcPfVLdlTW8eQ1ozltsP+eRBa6k2Obaae40eSS4OZ+NJMlH83goq6l3HPztb6RUV9b5250bU6TiGqmdbuq+e0/n+G77m8YP+FiBg0+dLO71k55kjDaIW40uVhJM6vj7703g1O/uIkTHG5s1TOQPcMhq7T9r2VRKt6ijJ2V26t4+MlneMLzAOniQWa/Cd0DaicaOzGjycUqmroRV0DgGGP4039W4Z47k7PS3NjxgidgVFS8xuQrFQ9Rxs4y+yCu/teX3CzLfIklVO1EYydmNLlYRbjqeFDguL//Fr9amMWL88r5+dBTsW16O/RZljYXqVQRRex4bWn83nsvmc6hTLzgMuStN8LXTjR2YkKTi1WEq44HBI7xuHjzzZd5cfsZ3DZuAD88ayCyeYCeZanUFkXseD2GE+wreHDSDfTMy4ZOWjtpa5pcrCJcddwfOMbjot7YeXlXH3478RiuHFN4aD0NDpXKIsSO156G121w42DCRZfRKy/70DoaN21Kk4uVhDrgC0rZeP5LzHrnFcoaBvOjay5n3KBu8SmfUlYVInY2ZR3Nr819jLItY8LE79H7mO/EqXCpSZOLVfk7Ipc4juGq2YZM53d5+pbjGNarU7xLppS1lc+jauXH/G5+RxaaYu6adA29e3SId6lSjiYXK/J3RHrd9RQbB6fl/Ia7br6Go3IzD32u7cVKHal8Hmb6BWS5XfzROCi/4CUGNiYWjZt2pcnFijaU4XXXY8OLU9z8fnQlGYGJJdmvIFaqhVxrPsXudmHHS7rNw8DaRcDpGjdxoHeitBhjDM/vKKTeOPBgw+ZIJ6M4oK041LBLpRQut5ffrcjHZRx4xX74pKUaN+1Oay4W4vUa/u/d5UxbkEXF0D9zS9E2pG9QFV6vIFbqCF6v4eevLWbGpm6cfNpTnJ65KuSoS42b9qPJxSKMMdzz1je8OK+cm07uy63nDUFEjlxQryBW6jDGGB6atYIZi7by87MHcfq4AUcupHHT7jS5WIAxhgffXcGL88r50bj+3HnWoNCJpZGO0VfqoKlz1zN17nquO7GIH57aP/yCGjftSvtcLOBPH6zmqc/Wc/1JRU0nFqXUQe8s3spvZq7gvGN6ct/5QzV2LCQmyUVExovIShFZIyJ3x2KbqeJfZev4y4eruaykQIMjBWnstNwX6/ZwxyuLKS3qwuTvjcBm09ixklYnFxGxA48D5wBDgStEZGhrt5sKZi/ddvCs67ffPUYTS4rR2Gm51TuqmPTMAgrzsphyzWgy0uzxLpIKEouaSymwxhizzhjjAl4CJsRgu0lt6ZYKfvryYo4tzGXy90Zg17OuVKSx0wI7K+u47un5pKfZmXb9ceRmOeNdJBVCLJLLUUB5wOvN/vcOIyKTRGSBiCzYtWtXDHabYMrnQdlkKJ/Hjso6bpw+n85ZaUz5fomedaUujZ2mBMQNwIF6NzdMn8++GhdPX3ccvTtnxbmAKpx2Gy1mjJkCTAEoKSkx7bVfSwi4OtjYnTya/SBVdQW8fuuJdO2QHu/SKYtL2dgJcS+j2z6ysXxrJVOvPY6jj9J59qwsFjWXLUBBwOve/vdUo8D7Srjr6bp7Pn+6bCRDenaMd8lUfGnsRBJ0L6MP33uDj1fu4sGLjmbcYJ0Z3OpikVzmA8Ui0ldEnMDlwNsx2G7y8F8d7MWOyzg4ekBfzt77vO/MLKjar1KKxk4kjVfVix23OPh4k5tpA8q4qtcO3+caO5bW6mYxY4xbRG4D3gfswFPGmGWtLlkyKShl84Uv8eprL5LRqRu3bPkzbHKBzQ4IeN06mV4K0thpgv+q+hVfzGL615X8X/pzpG1xw/RpMP5hmH23TkRpYTG5zsUYM8sYM9AY098Y81AstplMal0ebvzQxrNpl3DV8Bzk4AR6DTqZXorT2IlsgWcAExaP4ZguHtJwI42xsmKGxo7F6RX67eC3s1awckcVf7psJB0HjztY1ceeFvBcJ9NTKtDaXdXc9MwCeudmcsGF30MCY2XIBI0di9O5xdrYB8t38OwXG7l5bF++M7Ar0PXwCfRAJ9NTKsjOqjqufWoeDpsw7fpSOuZlHTnxZPehGjsWpsmlDe2squOu15cwtGdH7jx70KEPgifQ08BQ6qAD9W5unLaAPdUuXv7B8RTm+a9lCRU3GjuWpc1ibcTrNdz56hIO1Lv5yxUjSXfohZJKNaXB4+W2F75i2dYKHr/qWIb3zo13kVQLaXJpI9P+u4E5q3Zx7/lDGdCtQ7yLo5TlGWO4+/Vv+HjlLh6aeAynDe4e7yKpVtDk0gZWbKtk9uy3mdzzA67utT3exVEqIfx+9kpe/2ozfxhTxxX1r+r1KwlO+1xiqXweDWvn8P4Xe3nG8STp+z3IMy/qGHylmvDOu29h+3wWjw3ow4Rlf9XrV5KAJpdY8c+DZHfX8yMjOMQgxntoDL4GiFIhffifdzhj/s2cm+bGttXmixuNnYSnzWKxsqEMr7seG15s4kVsNh2Dr1QT3vtmG199+g5OcWPHi3i9IBo7yUBrLjGyK6+UDjhw4MbuSPdNT1G7R8fgKxXGp6t2cftLX3NJtzHYqt4+1BSmsZMUNLlEq3xe2Au2XG4vN39sI9vcx19PqKbLsNM1KJSCsHHz3zW7+cGzCyju1oG7J52J7D5GL4hMMppcohF0X4ngTsaHZi5nUfl+nrjqEroc0zOOBVXKQsLEzWdrdnPj9Pn06ZLNszeW0ikzTS+ITELa5xKNgPtKBE+S9/birUz/fCM3ntyXczSxKHVIiLj5bM1ubpjmSywv3DyGvBy9WV6y0uQSjYD7SgR2Mn65bg93vrqYkj6dufucwXEupFIWExQ3C21Hc8O0+fTN18SSCrRZLBr++0oEtgkv31rJTdMXUNgliyevKSHNrnlaqcMExM2nrsHcONPDoB4deOaGUk0sKUCTS7QC2oQ37anhmqfmkZPh4JkbSumc7Yxz4ZSyqIJSpm7syoP/Wc4J/bow5ZrRdMhIi3epVDvQ5NJM32731Vg8Xi8vTTqBXlXfwDc6ykWpYG6Pl9+99y1T565n/LAePHb5SDLSAiZwjTACUyW+5E4u0Ry84ZYJ8f6sb7ZxxyuL6ZDh4JkbxjCgfkXEUWRKJaRof/QjxE7t6k/4w7ddeXpTN647sYhfnT8Uu00OX1djJ6klb3KJ5uANt0zQ+/VXvsljqzrzxCdrGVWYyz+uHk23jhlQFmIUmQaISmTR/uhHiB3vtAtI87i4yzg46bSpnHHWsCPXDzUCU2MnqSRvL3SE4cNNLhPwvtft4slnn+GJT9ZyRWkBL0463pdYIOwoMqUSVjRxE2Y5r9ew4NO38bpdOPCSYfNwRubq0Otr7CS95K25NB68jWdWmXlQNhlXwYlszDya6no3WZkjKLalIV7w2tLYUJ2O+9UH2FCbwak4sBtDA3b2dh/Dq+ecwHFFXXzbDmwOCL71qlKJLEzcHHF8By23nw68/9hP+Xq3jWPS07Dj9t3zPjBpBDejaewkNTHGtPtOS0pKzIIFC1q+gWa0CZv1ZaypdtJn/oPYTQMu4+Aq1y/5ygwEYJSs4njbCvaaHO5Pe5Y03DTg4J9ZNzOoYwMDS8+h/+jTDt+3thUnHRFZaIwpiXc5mtIusdO4TGYezL47/LHuv8XE3K1ejl/5CGm4we7Efu7DSO3ew/ehcZOUIsVN4tVcojxIPV7Dq9t68NTCEs7Y/Tw/S2vAjpd08fDgyP3sHHmcf0nf377L/0nGYjeCF7t4+MmJeTD2jiP3r23FKlFF+wPfOOy+bHLYY93jNcza25vffzmKCytf4pQ036zGmAao3Xtk7GjcpJzESy5RHKTle2v42SuLmL9hH0N7dmT0qRdg/9I366rN7mTYiecxrKDb4dvNOhuWPQ4e15HV+UDBzQbaVqwSRXN/4EMc6xW1Dbwyv5zpn29g875aBvfowPjzL8H+4duRY0LjJuUkXnJp4iCdsWgL9765FAP86bIRXDTyKEQEBjfRvhttG7C2FatE1dwfeP+xXr3yE770DuHlj22Urf6Q2gYPpUVduOfcIZw1rIdviHHvGMWXShpJ0+ficnv59TvLeP7LTYwqzOWxy46lMC8r6vVVakvVPheP13DA5eZAvZt9BxrYXV3Prqp6yvfV8O22KlZsr2TjnhoAenbK4LrCnUzIXUeP4Wdq7Kgk63OBI6bn3lNdz63PfcW8DXv5wSn9+PnZg3CEm+tLOxZVivF6DSt3VPHluj0s3JTG9ooT2PXFAXZXv091vTvkOiLQNy+bYb068r2SAk4d1JWh7m+RZ67yxc7Xf9XYURElTnIJU9tYVL6fHz3/Fbur63nsspFcdOxRkbejHYsqRdQ1eHh35gy2Lf4PH9cN5CszkF6dMijMy+KY3rnk5zjplJlGTrqD7HQHuZlp5HdIJz8nnR4dM8h02g/fYNlcjR0VNWsll4AEUt1tFC/N20RRXjanZK3D+dzEw2obuzuP4JHZK3llYTndO2Tw6i0nMLx3btP70I5FlYwCYsf0Po63F29l5qwZ/Ln+fpzi5tZMJ3sufpXuw05p+T40dlQzWCe5BDRXee1p3Jv+IG/t8dVCfprxDrdRjx0vHreLD2a+xp3b91Hr8nDTyX35n9OL6ZiRFl1finYsqmQTEDvG7uTxwsk8ujyXBzqvIEN8w+sxDXTfOx8Ik1w0dlSMtSq5iMilwAPAEKDUGNPynsbDplwx9HV/xdPXXwTAks/34Vr/uv8CRzvP7yigtKgLvzh3CAPql8P8vzR9wVcgvaWqirO2ix0XNSs/5Ufjfsb3B1+NPPtq+JpGtBdLBtLYUVFqbc1lKfBd4J+tLknRWDy2NIzb4BYHl158Bb0G+a5FGTfoajwbi5GNc8nsO5ZnQl31KwLG63toe7CyvpjGDnYnXrcLl7HTYfA4bjlrkG8IfriahsaOamOtSi7GmBWA7yBurYJSKr/3Ol98NINTzppIr/4nHlZVt/cZA33GHL5OYOe8sYHNBoi2ByvLi3XszD1xKv/9aAbewpP4+Sn9kLl/PJRQQiUKjR3Vxtqtz0VEJgGTAAoLC0Mu03nQyZwz6GTfi2iGDAd3MI5/GGr3aHuwSirRxM7LO3qxo/f1PHsm2J+d0HQTl8aOamNNJhcR+QDoEeKje4wxM6LdkTFmCjAFfBeCNblCNEOGtYNRWVh7xs5jl42ktsFD+rw/RzdcWGNHtbEmk4sx5oz2KMgRoh32qB2MyqLaM3bsNiEn3dG84cIaO6oNWWcocjA9s1Kq+TRulEW0dijyROCvQFdgpogsMsacHZOSgZ5ZqaTVprGjcaMsIC4TV4rILmBjM1fLB3a3QXFaymrlAeuVyWrlgfBl6mOM6drehWmuFsROIv0fxIvVygPWK1Oz4yYuyaUlRGSBlWattVp5wHplslp5wJplaktW/L5WK5PVygPWK1NLyhNm6mCllFKq5TS5KKWUirlESi5T4l2AIFYrD1ivTFYrD1izTG3Jit/XamWyWnnAemVqdnkSps9FKaVU4kikmotSSqkEoclFKaVUzCVMchGRB0VkiYgsEpF/i0gvC5TpERH51l+uN0UkN87luVRElomIV0TiOoxRRMaLyEoRWSMid8e5LE+JyE4RWRrPcsSL1WLHanHjL5MlYsdKceMvT4tjJ2GSC/CIMWa4MWYk8C5wX5zLA/Af4GhjzHBgFfCLOJen8R4hc+JZCBGxA48D5wBDgStEZGgcizQNGB/H/ceb1WLHanEDFogdC8YNtCJ2Eia5GGMqA15mA3EfiWCM+bcxxu1/+QXQO87lWWGMWRnPMviVAmuMMeuMMS7gJWBCvApjjJkD7I3X/uPNarFjtbgBy8SOpeIGWhc71p24MgQReQi4BqgAxsW5OMFuAF6OdyEs4iigPOD1ZmBMmGVVO7Bw7GjcHJJUcWOp5NLU/S+MMfcA94jIL4DbgPvjXSb/MvcAbuB5K5RHpR6rxY7V4ibaMqnYsVRyacb9L54HZtEOyaWpMonIdcD5wOmmHS4aitv9dZpnC1AQ8Lq3/z3VRqwWO1aLm2jKZAFJFTcJ0+ciIsUBLycA38arLI1EZDxwF3ChMaYm3uWxkPlAsYj0FREncDnwdpzLlLKsFjsaN2ElVdwkzBX6IvI6MAjw4pty/BZjTFyzuoisAdKBPf63vjDG3BLH8gTeI2Q/ENv76zSvLOcCjwF24CljzEPxKIe/LC8Cp+KbNnwHcL8xZmq8ytPerBY7VosbsE7sWClu/OVpcewkTHJRSimVOBKmWUwppVTi0OSilFIq5jS5KKWUijlNLkoppWJOk4tSSqmY0+SilFIq5jS5KKWUijlNLkoppWJOk4tSSqmY0+SilFIq5jS5KKWUijlNLkoppWJOk0sciMh1IjI31ssqlWw0VhKXJheFiDhF5DUR2SAiRkROjXeZlLIiESnyx0h1wONXAZ+ni8hTIlIpIttF5GfxLG88WepOlCqu5uK7j8SrcS6HUokg1xjjDvH+A0Ax0AffLZU/FpHlxpjZ7Vk4K9CaSxsSkbtFZK2IVInIcv8NiUItZ0TkdhFZJyK7ReQREbEFLfOoiOwTkfUick7A+9eLyAr/PtaJyA+aW05jjMsY85gxZi7gafYXVaqVEiVWonAt8KAxZp8xZgXwJHBdG+zH8jS5tK21wFigE/Br4DkR6Rlm2YlACTAK361obwj4bAywEt/d4P4ATBUR8X+2E9+9yDsC1wN/EpFRACJSKCL7IzyujOm3VarlEi1WNorIZhF5WkTy/dvoDPQEFgcstxgY1vx/jiRgjNFHOz2ARfiC4TpgbsD7Bhgf8PqHwIf+59cBawI+y/Iv3yPMPt4CftyKMm4GTo33v5U+Uvth1VgBcvAlNgfQHXgNeN//WYF/fxkBy58JbIj3v2c8HlpzaUMico2ILGo8+wGOxndGFUp5wPONQK+A19sbnxhjavxPc/z7OEdEvhCRvf59nBthH0pZUqLEijGm2hizwBjjNsbsAG4DzhKRDkC1f7GOAat0BKqas49kocmljYhIH3ztrbcBecaYXGApIGFWKQh4XghsjWIf6cDrwKNAd/8+ZjXuw1/Vr47wuKqFX0+pmEnwWDH+vzZjzD5gGzAi4PMRwLKmypeMdLRY28nGd+DtAl9nIr6zsXB+LiJf4jvL+jHwxyj24QTS/ftw+zsvz8IXmBhjNvm31yR/8DUGs1NEMoB646/bK9WGEiZWRGQMsB9YDXQG/gJ8Yoyp8C/yDHCviCzA12x2M77+nZSjNZc2YoxZDkwGPgd2AMcAn0VYZQawEF9b80xgahT7qAJuB14B9gFXAm+3sMgrgVrgKOB9//M+LdyWUlFLsFjpB8zG19S1FKgHrgj4/H58gxM2Ap8Cj5gUHIYMIHpiGn8iYoBiY8yaeJdFKSvTWEkcWnNRSikVc5pclFJKxZw2iymllIo5rbkopZSKubgMRc7PzzdFRUXx2LVSIS1cuHC3MaZrvMvRFI0dZSWR4iYuyaWoqIgFCxbEY9dKhSQiG+Ndhmho7CgriRQ32iymQiufB2WTfX+VUtHT2AH0Cn0VSvk8mH4heFxgd8K1b0NBabxLpZT1aewcpDUXdaQNZb7gMB7f3w1l8S6RUolBY+cgTS7qSEVjfWddYvf9LRobflltAlDqkGhjJwXiRpvF1JEKSn3V+Q1lvuAIV63XJgClDhdN7KRI3GhyUaEVlDZ9wIdqAmhNkJTPazqhKWV1TcVOrOPGomLSLCYiT4nIThFZGovtqQQRqQmgudX+xrO5jx7y/U3i5oJGGjcpqqmms5bEjgWb2GJVc5kG/A3fvQxUqgjXBNCcan9jbaVic0qczQWZhsZN6onUdNbcJjMLN7HFJLkYY+aISFEstqUSTKgmgGir/YGBYbODzQFemh5EkCQ0blJYuKaz5sSOxU/K2q3PRUQmAZMACgsL22u3Khaa2xfSWO1vPJsKlygCA8kLjL4GOhVon4tKXdHEToKclLVbcjHGTAGmAJSUlOhUzImiJdXuaEebBQfSiCsPb1rTzn1AT8wSWnOP42hip6mTMovEjo4WU5G1dGRLNKPNYtFnkwL0xCxBtfQ4bip2mjops0jsaHJRkUXbxNVSremzUcrK2uo4jlS7sVDsxGoo8ovA58AgEdksIjfGYrvKAhoP5NPuab+zoOChmpl5lhxq2VoaN0muOTNdNFdBKYy948h4tFDsxOVOlCUlJUanDU9QkdpzY9nW27itzDyYfXebV/NFZKExpiTmG44xjZ0EE4/+j3aMnUhxo81iqS7w4IeWT1sRbVtvtMHW2FxWNvnwav7iFyzRWalUVAKbfcMd+7FOQOFip52byDS5pIJIB3XgkEYEvO7wySFSe240bb0t6WwM7POx2eHrFyKXUan2EBxTTSWIcMd+Sy44buklAY1NZO10YqbJJdlFOngPSwhe/womfHKI1LkfTcd/SzobAzsvKzbDwunNPxOzyNBMlYBCHTvBMTX+4aabn8Id+y254LgllwS0U/NyIE0uyS7SwRtcKwisuYRKDpFGqTQ1pUXjAd6SkWeN1fzyebDoxeatb6GhmSrBhDt2gmNqxYymE0S4k6+WXHDckksCWtpE1ooTM00uyS7UwRt4wAQmBGj6QIo0Bj/UZ6HO8mr3tKwWEe3FmYE2lGE8LsQCQzNVggn3gx4cU0MmwMbPIyeIcMduU8d0a0/MGrXkkoJWnphpckl2wQcvHHnAjL3j8OVjKThAa/ccvr9mMr2PozL/WPYecLFnw172HHCx74CLvTW+v/tqGthf08D+Ghf7axsoPJDB4147aWJwWGx6DGVxkfosghNC96EtPzEL974FTsxaMyBAk0sqCDx423sESZRnTB6vYVdVPdsr69hRWcfOyjp2VtWzq6r+4N/d1fXsqXbhOtg/dLiMNBuds5zkZjnJzUyjuFsOuVljedX9OCPcSxkx9nyttajoNdVnEXxSZvETs2aXsZUXUGtySSbRtI+29RX3wfwBWrf6U7Z2Hs3qykK2zF3P1v21bKuoY2tFLdsrfInE4z38miubQH5OOl07+B6DenQgL8dJfnY6eTlO8nLS6ZLlpEuOky5ZTjJ3LAzz/YcD32vb76mSUzyH9bZ3rAZrSW0ngCaXZBFt+2grD5hwvF7D9so6Nuw+wIY9NWzcc4CNe2oo31fDpr01VNUNBWqBhYCvltErN5NenTI5aUA+PTtl0L1jxsG/3Tqmk5edjt0mLfv+rWlCUCpYPH7o2yhWjxDppLQVNTJNLsmiOe2jrThg6ho8rN1VzZqd1azdWc3aXQdYu6uaDXsOUNdwqLnKabfRu0smhV2yGN2nMwWds+jdOZPenbM4qnMmnbPSEIkycUQS6r4W7nqYdQcYoyPEVHSaqvW31w99qP225b7a8KRMk0uyiPGZlddr2Li3hm+3VbJiWyXfbq9i9c5qNu45QGPrlU2gsEsW/bvmcPKAfIrys+mXn02f/Gx6dMyIvtbRUuHuayECxut76Agx1ZTm1PqT5Thqh5MyTS7JohVnVh6vYe2uahaX72fZ1kqWbqlgxbZKDrg8gC+JFOVnM7hHBy4Y0YuB3XMo7taBovws0h32tvpGTQt3X4vgzlcdIaYisdBMwu2inU7KNLkkk2jmMQJ2VtXx1cb9fF2+j6837Wfplgpq/Ikky2lnaM+OXDK6N0N7dWRIz44Ud+tApjOOSSScSPe1iGZoqFIQ/47z9tZOJ2WaXBJdE9NTGLuTbRNeZk5tX+Zv2Mf8DXvZtLcGgDS7MLRXJy4d3ZvhvXMZUdCJfvk52Nq6OStWmpoxQJOKikaoa8HacQ6udtdOJ2WaXKykuVMthGkrrljxMR3c9djw4nHX89xLz/N3zwTysp2UFHXmmhP6MKpPZ4b27EhGmgVrJM2hSUTFQuAUQ8k+XVA7nZRpcrGKlhzUAdVbr8fFh7Ne47dVNXTe4+R5p4M0ceOVNEaOPZ8PR3+HfvnZsRmhpVSiC3cilyr9L+1wUqbJxSoiHdQhAmFbRS2Lqos5DQd2Y2gwdp7a3JuiflmMPf4idnU4loKKhTj6juWsZAwOpQI1p9Yf6UQu1fpf2pAmF6sId1AHBILXnsabxzzB9PJuLNlcAdg5N/fXXNxlA3lHn8bTx50R0MzVFzg1Pt9FqfbU3Fp/pBO5eF3PkoQ0uVhFmIO6YsVHB/tPvA2GtfNmI72u5a7xgzhraHf6d83Rpi6V2ppZ62+ydqL9eDGhycVK/Ad1VV0D787bxBtfbcazMf1g/wn2NK6/7PvcNfSkeJdUKeuIotZ/WI1GayftQpOLRRhj+GLdXl5ZUM57S7dR1+Clf9dsJp55PhXdRtNj33woGktXDQSlDhcuWTTV/KWx1KY0ucRZRU0Dr321mee/3Mi6XQfokOHg4lG9uaZgBwNrFiJ9s6DgFOCUQyvpbXuVOlyoZKGd83GlySVO1uys4unPNvD6V5upa/BybGEuj146gvOO6embOn761b6gmBPUQZkK4/CVigVt/oorTS7tyBjD5+v28I9P1zFn1S6cDhsXjezFtScWMaxXp0MLRqrOp8o4fKVaKrhmH8WUSCr2NLm0A6/X8J8VO/j7J2tZXL6f/Jx07jhzIFeOKSQvJ/3IFSJV57Wqr1R44Wr2Qe97vz8DV68SXB4vLreXBo+XBrfB5fE9d3t8z90eL25v43ODx+ulwWNwextfG9zeQ38bl/f6X3uNOey1x/+e57DnHHrPGIxp/Mz32+E1/ufG/9y/vPG/5/F/bvyf+973vTYB6zU+N/gmPjb4tgWE/Bx82/352YO4orSw2f8VmlzakDGG95ft4LEPVvHt9ioKu2Tx0MSjuXhU78jTrjQ1PYNW9VUSM8ZQ1+Clqq6Bqno3VXVuaurdVNe7OeByU+PyUFPv4YDLTW2Dh1qXhxqXh7oGD+N2vchF7nrseHG765n+3DNMtx/g8vpXmeSpx4EXd0M9f5wylb97drf5dxEBh02w2wS7CDb/c4dNEPG9Z7cJNhsHP7eJBDwHu39Zm3DwMxGw2SDNZsMmjdsCEUEAm833VwLWF/Avi/8z/3aCPxMA3/5EoE9eVou+uyaXNmCM4ZNVu3j0/ZUs21pJ3/xs/nTZCC4Y3guH3RbdRiKNZtGRLiqB1Ls97K52sbOyjl1V9ew54GLvARd7ql3sr3Gxr8bFvpoGKmsbqKhtoLKugQaPaXrDQLrDRqbTTmaa7+GUIZyPA3DjkTT25JcyqkMuxnUy3vWv4TENGHsafUacxV1dBuG023A6bKTZGx9Cmt2G027D4X/usAkO/2cOm//vwfcPJY7G9xqTR+OPeqrS5BJji8v387v3VvDFur0Udsli8qUjmDCyGUlFqQRT1+Bho//W1pv21rB5Xy2b99WyvbKWbfvr2HPAFXK9nHQHuVlpdM5ykpuVRu/OmXTMTKNTZhodMhx0yEijY4aDbKeDnAwHOekOstMdZDvtZKU7yEyzh7gh3alQPhI2lGEvGstdB0/CjoXyAQffv0xPztpcTJKLiIwH/gzYgX8ZYx6OxXYTyfaKOn733gpmLNpKl2wnD1wwlCvH9MHp0KSikoPHa1izs5pvt1eyYlsV326vZM3Oarbsr/W30fvkpDs4KjeTnrkZHHNULj07ZdCtQzpd/Y/8nHS6ZDvbbkbucDV7rfG3q1YnFxGxA48DZwKbgfki8rYxZnlrt50I6ho8TJ27nsc/XoPba/jRuP7c8p3+dMhIi3fRlGqVqroG5q3fy7wNe1m0aT/fBNxULs0u9O+aw7GFnbl4VG/6dc2mb342BZ2zyM1KS+nmIOUTi5pLKbDGGLMOQEReAiYASZ9c/rtmN/e8tZT1uw9w9rDu3HveUAq6hOn80mGQKgQr1fqNMSzdUsl/lm/n09W7WbqlAo/X4LTbGNKrI5eO7s2IglyG9epEv67ZpGlTr4ogFsnlKKA84PVmYEzwQiIyCZgEUFjY/GFtVrL3gIuHZq7g9a820ycvi2dvLGVscdfwK+iFjyoEK9T6jTEs2VzBG19tZvay7eyorMcmcGxhZ354an9O6J/HqMLOiX9TOdXu2q1D3xgzBZgCUFJSEt1QEAt675tt3PvWUipqG/jRuP78z2nFTQeeXvioQotbrX9/jYsX55Xz2sJy1u46gNNh47RB3ThzaHfGDe5Gl2xnWxdBJblYJJctQEHA697+95LKvgMu7n97GW8v3srRR3Xk+ZvHMLhHx+hWjnThozaXpbKoav2xtG5XNU99tp7XFvqmHTquqDM3j+3HOcf0pFOmBfsJNT4SViySy3ygWET64ksqlwNXxmC7ljF39W5++soi9h1w8dMzBvLDcf2b194cfOEjQNlkyMyD2Xdrc5mKKBZNylv21/LHf6/ija83k2b3TTt0w8l9oz9BigdtTk5orU4uxhi3iNwGvI+vU/IpY8yyVpfMAho8Xib/exX/nLOW/l1zePq64zj6qE5NrxhK4zDIwIARAeP1PbS5LBVFVetvTZNyZV0Df/1wNdM/3wjAzWP7cfPYfnTtEGLaIavR5uSEFpM+F2PMLGBWLLZlFZv31XDbC1+zqHw/V5QWct/5Q8l0xqBTMzBgjM03hwOi84Slpjat9X/07Q5++cZSdlTVcfGo3vzszIH0ys2M1ebbns6jl9D0Cv0QPl21ix+/9DUer+HxK0dx3vCesdt4cMCMfxhq92ibcgpqq1r//hoXv35nOW9+vYWB3XP45/dPYkRBbms32/7CNSc3xor2x1iaJpcAXq/h8Y/X8McPVjGoewf+cfVoivKzY7sTnXhSBYh1rX/51kp+8NwCtu2v4/bTi/nRuP6kOxJ4GHGo5uTGkzLtr7Q0TS5+NS43P3t5MbOXbeeikb347XePIcvZRv88Og2FagMzFm3hf19fQm6mk1dvOYFjCzvHu0ixE9z/smKG9sdYnCYXYFtFLTdNX8CKbZXce94Qbjy5r05foRLK72at4J9z1lFa1IXHrxqVGB32zRHcnDxkAmz8XPtjLCzlk8vi8v3c9MwCal0epl57HOMGd2v5xgLbgEGbvlS76ZiZxnUnFnHPeUOSc1qWUM3J3YdqjFlYSieXj7/dyQ+f/4q8HCfP3TiGQT06tHxjgW3CNjsg4HVre7BqFz88tX/y17aDm5O1ednSkvAUJzqvLijnpmcW0K9rNm/88MTWJRYIahNuOLI9WKk2lPSJRSWclKu5GGP4+ydreeT9lYwtzueJq0eTkx6Df4bANuHgmou2BysVPR1inBRSKrkYY/jD+yt54pO1TBjZi0cuGRG7m3mFGpOvAaJUdBoTik6JlDRSJrl4vYb/e3c50/67gavGFPLghKOxHXGL1AiiOZsK1SaslIpMp0RKSimRXDxewy/eWMIrCzZz08l9uee8Ic1ro9YJ9JRqOzolUlJK+uTi9Rrufn0Jry7czO2nDeCnZw5sfuenTqCnVNvRKZGSUlInF6/X8Ms3v/ElltOL+dmZA1u2IZ1AT6m2o1MiJaWkTS7GGH41YykvzS/nR+P689Mzilu+MT34lWqZaEd+6TUrSScpk4sxhodmruD5Lzdxy3f6c+dZg1p/HUC4g1+HTSoVmvZVprSkTC5/+2gN/5q7nnuHV3Fj9kxkc2X4gzpccogmaWjwKBWe9lWmtKRLLtM+W8/k/6zip4P2ceO6O5HVLpgT5oc/XHKINmlo8CgVnvZVprSkmv5lxqItPPDOcs4c2p3b+m1HmpqCJVRyiPR+sMbgEbsGj1LBGvsqT7tHa/UpKGlqLnNX7+bOVxdzfL8u/PWKY7Fv90LZo4fOmjLzDr+LHRx5ZtW4TGZe5DOuwCYz7ehXKrxIHfXaX5nUEjO5BB2US7dU8INnF9C/aw7//H4JGWn2w0d4hZtSItIy4cbah2oyG3tH/P4tlIq1aH/0W5MctL8y6SVecgk6KHdMfIXr33LRKTONadeX0ikz7dCyjQmkbHL4vpFwy9TuCZ00tJ9FJbNof/QjLRdN0tE4SnqJ1+cScFAaj4t3ZrxCfYOH6TeU0qNTRuh1oukbibb/RPtZVDKLtr8x3HKNSeejh3x/y+eFXl/jKOklXs3Ff1AajwsXDt6vGcCUG0oo7h7hfiyhZiwO7n+J9kJJvaBSJbNw/ZDBx3q4kWCRaiTBNRqNo6SWeMmloBRzzQzee/c1/lXei6svvYTj++VFtV6Tw4yjvUpYryZWySqavsrg5QJP2MINhgkXdxpHSStxkkvAWc/f1+bxyKZT+fHpxXx3VO/mbUfbelUqa86tIyL1VQYuF5w4Qg2G0bhLOYmRXAIOXo8tjQ9r7mbiseP4SUvmC9MLu1Sqau4IrWhjJThxhBoMo3GXcqyVXMKdVQV24rsN383bwKUXH3PkfGHRnpVpW69KRc2tPUQbK9EkDo27lNOq5CIilwIPAEOAUmPMghZvLNJZVdFYvPY0vG6DWxycf+GlpDvs0a8fTNt6VSpqSe0hOFZCncA1ZzCMxl3KaG3NZSnwXeCfrS5JhLOqmu6juC/rQXpXfMXEiZfTZ+DJh9ZrPNgrNmubrlKRtLb2EIvBMCpltCq5GGNWAK2fzh5Cn1WVz8O7voy/rezKG7uOYuq1F9FncLdD6wQe7DY72BzgRdt0lWoUXNNoTRLQTnnVDO3W5yIik4BJAIWFhUcuEGpo4/QLwV3P/xgHg0+awrjAxAKHH+xeYPQ10KlA23SVgthPsaKd8qoZmkwuIvIB0CPER/cYY2ZEuyNjzBRgCkBJSYkJuVDgWVXZZLzuemx4cYqbCzquO3L54IN9xJWaVFRqCtUXEuuahnbKq2ZoMrkYY85oj4IEW5k5gkLjwClubI50pK+OQFEqpHA1lObUNPR2xCrGrDUU2W/zvhqumm04LuP/eLS0iuyBp+oIFJVUYjrSMlwNpamTr8aEEulKfKVaqLVDkScCfwW6AjNFZJEx5uzWbLOqroGbpi+g3u3ljklXk90twpxhSiWu2I20bGo+sKZmNRYB4/U9tKNexUhrR4u9CbwZo7Lg8Rpuf/FrVu+sZvr1pQzQxKKSVExHWkY7H1igwNqOsYHNBoh21KuYsVSz2OsLN/Pxyl385qKjObk4P97FUcoSmhxpCdHPB9YouLYT7uZ4SrWQpZLLJaN7k5fj5PQh3eNdFKVarV1HWjaKthNfB8OoNmap5GKziSYWlTTiMtKyOUlDB8OoNmSp5KKUigFNGsoCEu82x0olARGZKCKbgRPwjbR8P95lUiqWxJjITbhtslORXcDGZq6WD+xug+K0lNXKA9Yrk9XKA+HL1McY07W9C9NcLYidRPo/iBerlQesV6Zmx01ckktLiMgCY0xJvMvRyGrlAeuVyWrlAWuWqS1Z8ftarUxWKw9Yr0wtKY82iymllIo5TS5KKaViLpGSy5R4FyCI1coD1iuT1coD1ixTW7Li97VamaxWHrBemZpdnoTpc1FKKZU4EqnmopRSKkFoclFKKRVzCZNcRORBEVkiIotE5N8i0ssCZXpERL71l+tNEcmNc3kuFZFlIuIVkbgOYxSR8SKyUkTWiMjdcS7LUyKyU0SWxrMc8WK12LFa3PjLZInYsVLc+MvT4thJmOQCPGKMGW6MGQm8C9wX5/IA/Ac42hgzHFgF/CLO5Wm8R8iceBZCROzA48A5wFDgChEZGsciTQPGx3H/8Wa12LFa3IAFYseCcQOtiJ2ESS7GmMqAl9lA3EciGGP+bYxx+19+AfSOc3lWGGNWxrMMfqXAGmPMOmOMC3gJmBCvwhhj5gB747X/eLNa7FgtbsAysWOpuIHWxU5CTVwpIg8B1wAVwLg4FyfYDcDL8S6ERRwFlAe83gyMiVNZFJaOHY2bQ5IqbiyVXJq6/4Ux5h7gHhH5BXAbcH+8y+Rf5h7ADTxvhfKo1GO12LFa3ERbJhU7lkouzbj/xfPALNohuTRVJhG5DjgfON20w0VDcblHSPNtAQoCXvf2v6faiNVix2pxE02ZLCCp4iZh+lxEpDjg5QTg23iVpZGIjAfuAi40xtTEuzwWMh8oFpG+IuIELgfejnOZUpbVYkfjJqykipuEuUJfRF4HBgFefFOO32KMiWtWF5E1QDqwx//WF8aYW+JYnonAX4GuwH5gkTHm7DiV5VzgMcAOPGWMeSge5fCX5UXgVHzThu8A7jfGTI1Xedqb1WLHanED1okdK8WNvzwtjp2ESS5KKaUSR8I0iymllEocmlyUUkrFnCYXpZRSMafJRSmlVMxpclFKKRVzmlyUUkrFnCYXpZRSMff/vm68/9xyjf8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ind = ['alpha_%.2g' % alpha_ridge[i] for i in range(0, len(alpha_ridge))]\n",
    "coef_matrix_ridge = pd.DataFrame(index=ind, columns=col)\n",
    "\n",
    "models_to_plot = {1e-15: 221, 1e-3: 222, 1: 223, 50: 224}\n",
    "for i in range(len(alpha_ridge)):\n",
    "    coef_matrix_ridge.iloc[i, ] = ridge_regression(data, predictors,\n",
    "                                                   alpha_ridge[i],\n",
    "                                                   models_to_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rss</th>\n",
       "      <th>intercept</th>\n",
       "      <th>coef_x_1</th>\n",
       "      <th>coef_x_2</th>\n",
       "      <th>coef_x_3</th>\n",
       "      <th>coef_x_4</th>\n",
       "      <th>coef_x_5</th>\n",
       "      <th>coef_x_6</th>\n",
       "      <th>coef_x_7</th>\n",
       "      <th>coef_x_8</th>\n",
       "      <th>coef_x_9</th>\n",
       "      <th>coef_x_10</th>\n",
       "      <th>coef_x_11</th>\n",
       "      <th>coef_x_12</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha_1e-15</th>\n",
       "      <td>0.62</td>\n",
       "      <td>0.9</td>\n",
       "      <td>-0.12</td>\n",
       "      <td>0.56</td>\n",
       "      <td>1.2</td>\n",
       "      <td>-2.9</td>\n",
       "      <td>-3.8</td>\n",
       "      <td>1.4</td>\n",
       "      <td>3.7</td>\n",
       "      <td>1.2</td>\n",
       "      <td>-0.53</td>\n",
       "      <td>-0.46</td>\n",
       "      <td>-0.11</td>\n",
       "      <td>-0.0096</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_1e-10</th>\n",
       "      <td>0.66</td>\n",
       "      <td>0.94</td>\n",
       "      <td>-0.014</td>\n",
       "      <td>-0.095</td>\n",
       "      <td>0.24</td>\n",
       "      <td>-1.1</td>\n",
       "      <td>-1.3</td>\n",
       "      <td>0.31</td>\n",
       "      <td>1.1</td>\n",
       "      <td>0.65</td>\n",
       "      <td>0.1</td>\n",
       "      <td>-0.024</td>\n",
       "      <td>-0.01</td>\n",
       "      <td>-0.001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_1e-08</th>\n",
       "      <td>0.71</td>\n",
       "      <td>0.96</td>\n",
       "      <td>0.018</td>\n",
       "      <td>-0.58</td>\n",
       "      <td>-0.29</td>\n",
       "      <td>-0.021</td>\n",
       "      <td>0.29</td>\n",
       "      <td>0.24</td>\n",
       "      <td>0.08</td>\n",
       "      <td>0.013</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>-0.00022</td>\n",
       "      <td>0.00011</td>\n",
       "      <td>3.4e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.0001</th>\n",
       "      <td>0.74</td>\n",
       "      <td>0.95</td>\n",
       "      <td>-0.047</td>\n",
       "      <td>-0.45</td>\n",
       "      <td>0.079</td>\n",
       "      <td>0.034</td>\n",
       "      <td>-0.0081</td>\n",
       "      <td>0.0024</td>\n",
       "      <td>-0.00012</td>\n",
       "      <td>-0.00012</td>\n",
       "      <td>7.1e-05</td>\n",
       "      <td>-2.1e-05</td>\n",
       "      <td>2e-06</td>\n",
       "      <td>2.3e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.001</th>\n",
       "      <td>0.78</td>\n",
       "      <td>0.94</td>\n",
       "      <td>0.011</td>\n",
       "      <td>-0.42</td>\n",
       "      <td>0.018</td>\n",
       "      <td>0.01</td>\n",
       "      <td>-0.0047</td>\n",
       "      <td>0.0014</td>\n",
       "      <td>-0.00026</td>\n",
       "      <td>2.3e-05</td>\n",
       "      <td>9.2e-06</td>\n",
       "      <td>-6.1e-06</td>\n",
       "      <td>2.2e-06</td>\n",
       "      <td>-4.9e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.01</th>\n",
       "      <td>0.99</td>\n",
       "      <td>0.87</td>\n",
       "      <td>0.091</td>\n",
       "      <td>-0.29</td>\n",
       "      <td>0.018</td>\n",
       "      <td>0.00085</td>\n",
       "      <td>-0.0022</td>\n",
       "      <td>0.00084</td>\n",
       "      <td>-0.00025</td>\n",
       "      <td>5.8e-05</td>\n",
       "      <td>-8.8e-06</td>\n",
       "      <td>-6.7e-07</td>\n",
       "      <td>1.4e-06</td>\n",
       "      <td>-7.9e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_0.1</th>\n",
       "      <td>1.9</td>\n",
       "      <td>0.77</td>\n",
       "      <td>0.2</td>\n",
       "      <td>-0.12</td>\n",
       "      <td>0.018</td>\n",
       "      <td>-0.0033</td>\n",
       "      <td>0.00023</td>\n",
       "      <td>7.7e-05</td>\n",
       "      <td>-5.4e-05</td>\n",
       "      <td>2.1e-05</td>\n",
       "      <td>-7.1e-06</td>\n",
       "      <td>2.1e-06</td>\n",
       "      <td>-5.7e-07</td>\n",
       "      <td>1.4e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_1</th>\n",
       "      <td>4.6</td>\n",
       "      <td>0.59</td>\n",
       "      <td>0.13</td>\n",
       "      <td>-0.049</td>\n",
       "      <td>0.012</td>\n",
       "      <td>-0.0028</td>\n",
       "      <td>0.00064</td>\n",
       "      <td>-0.00013</td>\n",
       "      <td>2.4e-05</td>\n",
       "      <td>-2.8e-06</td>\n",
       "      <td>-3.8e-07</td>\n",
       "      <td>4.6e-07</td>\n",
       "      <td>-2.4e-07</td>\n",
       "      <td>1e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_5</th>\n",
       "      <td>9.8</td>\n",
       "      <td>0.41</td>\n",
       "      <td>0.056</td>\n",
       "      <td>-0.021</td>\n",
       "      <td>0.0058</td>\n",
       "      <td>-0.0016</td>\n",
       "      <td>0.00046</td>\n",
       "      <td>-0.00013</td>\n",
       "      <td>3.6e-05</td>\n",
       "      <td>-1e-05</td>\n",
       "      <td>2.8e-06</td>\n",
       "      <td>-8e-07</td>\n",
       "      <td>2.3e-07</td>\n",
       "      <td>-6.4e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_10</th>\n",
       "      <td>13</td>\n",
       "      <td>0.34</td>\n",
       "      <td>0.035</td>\n",
       "      <td>-0.014</td>\n",
       "      <td>0.004</td>\n",
       "      <td>-0.0012</td>\n",
       "      <td>0.00035</td>\n",
       "      <td>-0.0001</td>\n",
       "      <td>3e-05</td>\n",
       "      <td>-9e-06</td>\n",
       "      <td>2.7e-06</td>\n",
       "      <td>-8e-07</td>\n",
       "      <td>2.4e-07</td>\n",
       "      <td>-7.3e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_20</th>\n",
       "      <td>16</td>\n",
       "      <td>0.28</td>\n",
       "      <td>0.021</td>\n",
       "      <td>-0.0083</td>\n",
       "      <td>0.0025</td>\n",
       "      <td>-0.00076</td>\n",
       "      <td>0.00023</td>\n",
       "      <td>-7e-05</td>\n",
       "      <td>2.1e-05</td>\n",
       "      <td>-6.6e-06</td>\n",
       "      <td>2e-06</td>\n",
       "      <td>-6.2e-07</td>\n",
       "      <td>1.9e-07</td>\n",
       "      <td>-5.9e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_50</th>\n",
       "      <td>20</td>\n",
       "      <td>0.22</td>\n",
       "      <td>0.0095</td>\n",
       "      <td>-0.0039</td>\n",
       "      <td>0.0012</td>\n",
       "      <td>-0.00038</td>\n",
       "      <td>0.00012</td>\n",
       "      <td>-3.6e-05</td>\n",
       "      <td>1.1e-05</td>\n",
       "      <td>-3.5e-06</td>\n",
       "      <td>1.1e-06</td>\n",
       "      <td>-3.4e-07</td>\n",
       "      <td>1.1e-07</td>\n",
       "      <td>-3.3e-08</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rss intercept coef_x_1 coef_x_2 coef_x_3 coef_x_4 coef_x_5  \\\n",
       "alpha_1e-15  0.62       0.9    -0.12     0.56      1.2     -2.9     -3.8   \n",
       "alpha_1e-10  0.66      0.94   -0.014   -0.095     0.24     -1.1     -1.3   \n",
       "alpha_1e-08  0.71      0.96    0.018    -0.58    -0.29   -0.021     0.29   \n",
       "alpha_0.0001 0.74      0.95   -0.047    -0.45    0.079    0.034  -0.0081   \n",
       "alpha_0.001  0.78      0.94    0.011    -0.42    0.018     0.01  -0.0047   \n",
       "alpha_0.01   0.99      0.87    0.091    -0.29    0.018  0.00085  -0.0022   \n",
       "alpha_0.1     1.9      0.77      0.2    -0.12    0.018  -0.0033  0.00023   \n",
       "alpha_1       4.6      0.59     0.13   -0.049    0.012  -0.0028  0.00064   \n",
       "alpha_5       9.8      0.41    0.056   -0.021   0.0058  -0.0016  0.00046   \n",
       "alpha_10       13      0.34    0.035   -0.014    0.004  -0.0012  0.00035   \n",
       "alpha_20       16      0.28    0.021  -0.0083   0.0025 -0.00076  0.00023   \n",
       "alpha_50       20      0.22   0.0095  -0.0039   0.0012 -0.00038  0.00012   \n",
       "\n",
       "             coef_x_6 coef_x_7 coef_x_8 coef_x_9 coef_x_10 coef_x_11 coef_x_12  \n",
       "alpha_1e-15       1.4      3.7      1.2    -0.53     -0.46     -0.11   -0.0096  \n",
       "alpha_1e-10      0.31      1.1     0.65      0.1    -0.024     -0.01    -0.001  \n",
       "alpha_1e-08      0.24     0.08    0.013   0.0004  -0.00022   0.00011   3.4e-05  \n",
       "alpha_0.0001   0.0024 -0.00012 -0.00012  7.1e-05  -2.1e-05     2e-06   2.3e-06  \n",
       "alpha_0.001    0.0014 -0.00026  2.3e-05  9.2e-06  -6.1e-06   2.2e-06  -4.9e-07  \n",
       "alpha_0.01    0.00084 -0.00025  5.8e-05 -8.8e-06  -6.7e-07   1.4e-06  -7.9e-07  \n",
       "alpha_0.1     7.7e-05 -5.4e-05  2.1e-05 -7.1e-06   2.1e-06  -5.7e-07   1.4e-07  \n",
       "alpha_1      -0.00013  2.4e-05 -2.8e-06 -3.8e-07   4.6e-07  -2.4e-07     1e-07  \n",
       "alpha_5      -0.00013  3.6e-05   -1e-05  2.8e-06    -8e-07   2.3e-07  -6.4e-08  \n",
       "alpha_10      -0.0001    3e-05   -9e-06  2.7e-06    -8e-07   2.4e-07  -7.3e-08  \n",
       "alpha_20       -7e-05  2.1e-05 -6.6e-06    2e-06  -6.2e-07   1.9e-07  -5.9e-08  \n",
       "alpha_50     -3.6e-05  1.1e-05 -3.5e-06  1.1e-06  -3.4e-07   1.1e-07  -3.3e-08  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.options.display.float_format = '{:,.2g}'.format\n",
    "coef_matrix_ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
